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ABSTRACT 


Effective  guidance  laws  that  are  optimal  for  tactical  air-to-air  scenarios  tend  to  improve 
the  perfonnance  characteristics  of  the  missile  and  increase  the  probability  of  a  hit  in 
combat.  Proportional  guidance  is  the  current  baseline  algorithm  for  tactical  missile 
guidance.  Increases  in  computational  capabilities  now  permit  more  complicated  guidance 
laws  to  be  implemented.  This  research  focuses  on  two  promising  advanced  guidance 
laws,  comparing  them  to  proportional  navigation  using  simulation,  with  the  kinematic 
boundary  as  the  performance  measure.  Studies  are  also  made  of  perfonnance  degradation 
in  the  presence  of  sensor  noise. 

The  three  guidance  laws,  Proportional  Navigation  (PN),  Augmented  Proportional 
Navigation  (APN)  and  Differential  Geometry  (DG),  were  each  simulated  against  a  non¬ 
maneuvering  target  and  a  maneuvering  target.  The  theoretical  missile  engagement 
envelope  (the  kinematic  boundary)  is  utilized  as  a  simple  and  intuitive  visual  aid  in 
comparing  the  effectiveness  of  each  guidance  law. 

Band-limited  white  noise  is  then  introduced  into  the  seeker  system  to  detennine 
the  ability  of  the  guidance  law  to  deal  with  noise  perturbations,  in  particular,  to  discover 
the  level  of  noise  tolerance  for  each  guidance  law. 

This  research  used  a  simulation  model  previously  developed  here  at  the  Naval 
Postgraduate  School  (NPS).  This  simplified  six  degree  of  freedom  (6DOF)  model  was 
used  in  a  slightly  modified  fonn  to:  1)  verify  earlier  results  obtained  at  NPS,  2) 
investigate  an  additional  guidance  law,  the  DG  law,  and  3)  study  the  effects  of  noise  on 
the  robustness  of  the  various  guidance  laws. 
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I.  INTRODUCTION 


The  history  of  manned  flight  is  a  relatively  short  one.  In  the  span  of  one  hundred 
years,  we  have  seen  technological  innovations  and  courageous  pioneers  who  have  pushed 
the  development  of  aeronautics  to  where  we  stand  today.  From  the  day  the  Wright 
brothers  achieved  powered  flight  in  1903,  to  the  day  Major  Charles  E.  Yeager  broke  the 
sound  barrier  in  1947  [1],  to  the  unmanned  drones  that  fly  over  Afghanistan  today, 
aviation  history  is  full  of  stories  of  human  courage  and  ingenuity:  Courage  to  step  into  the 
unknown,  risking  life  and  limb  for  the  advancement  of  science.  Ingenuity  to  design  and 
develop  the  technologies  required  for  the  wonders  of  modem  day  flight:  planes  that 
takeoff  vertically,  flying  wings  that  look  more  like  a  child’s  boomerang  than  a  menacing 
deep  penetration  bomber,  missiles  that  can  be  launched  at  enemy  targets  even  before  the 
human  eye  can  see  them. 

It  was  perhaps  inevitable  that  with  the  advancement  in  fighter  aircraft  technology, 
something  more  than  simple  .50  caliber  rounds  that  fire  at  a  fixed  point  ahead  of  the 
aircraft  was  needed  to  shoot  them  down.  Rockets  have  been  known  for  hundreds  of  years, 
far  longer  than  human  flight,  yet  it  was  not  until  recent  decades  that  developments  in 
technology  allowed  for  the  evolution  of  the  simple  rocket  into  the  sophisticated  missiles 
of  today. 

In  this  thesis,  the  focus  is  on  tactical  missiles.  Tactical  missiles  are  used  in 
scenarios  where  the  ranges  concerned  are  more  limited  and  are  usually  guided  by  a  seeker 
sensor.  The  seekers  can  be  active,  passive  or  even  semi-active,  utilizing  electromagnetic 
waves  from  a  radar,  an  infra-red  (IR)  sensor  or  a  laser.  This  suite  of  sensing  abilities 
allows  the  seeker  to  detect  and  identify  the  target  and  guide  the  missile  to  it.  The  ability 
to  guide  a  missile  to  a  detected  target  is  the  province  of  guidance  laws.  Perhaps  the  most 
intuitive  and  also  one  of  the  earliest  guidance  laws  is  the  pursuit  guidance  law.  Pursuit 
guidance  basically  states  that  as  long  as  the  missile  is  pointed  at  the  target  at  all  times, 
given  enough  kinetic  energy,  the  missile  will  hit  the  target.  While  a  simple  guidance  law 
to  implement,  in  reality,  it  does  not  work  so  well  because  the  kinetic  energy  available  to 
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such  missiles  is  limited.  The  rocket  propellant  boosts  the  missile  up  to  maximum  speed 
within  seconds  of  launching  before  being  consumed,  and  the  missile  than  glides  the  rest 
of  the  distance  on  kinetic  energy  alone.  Hence,  in  pursuit  guidance,  the  missile  is  more 
apt  to  run  out  of  kinetic  energy  before  it  can  successfully  close  with  the  target.  This  is 
especially  true  when  launched  at  a  maneuvering  target  from  its  frontal  hemisphere. 

The  solution  to  this  problem  is  the  proportional  navigation  (PN)  guidance  law. 
While  not  as  intuitive  as  the  pursuit  guidance  law,  it  is  still  a  simple  and  robust  concept. 
Basically,  the  concept  of  proportional  navigation  is  to  point  the  missile  at  a  point  ahead  of 
the  target  so  that  the  missile  will  lead  the  target,  and  this  will  reduce  the  amount  of 
maneuvers  necessary,  thus  conserving  kinetic  energy  for  the  missile  to  make  the 
intercept.  The  implementation  is  simple  and  the  basic  idea  is  to  strive  to  maintain  a 
constant  line  of  sight  (LOS)  angle  between  the  missile  and  the  target,  the  premise  being 
that  a  constant  LOS  angle  would  signify  that  the  missile  and  the  target  are  on  a  collision 
course.  The  acceleration  commands  are  theoretically  applied  perpendicular  to  the  LOS 
and  are  proportional  to  the  LOS  rate  and  closing  velocity.  This  basic  concept  is  so 
successful  that  most  of  the  more  successful  guidance  laws  in  use  today  tend  to  be 
extensions  of  the  basic  PN  law. 

In  particular,  the  Augmented  Proportional  Navigation  (APN)  law  [2],  [3],  and  the 
Differential  Geometry  (DG)  law  [4],  [5],  are  investigated  in  greater  detail  in  this  study  as 
extensions  of  the  basic  PN  law.  Guidance  laws  based  on  Optimal  Control  Theory  [6]  are 
not  investigated  in  this  study  due  to  time  and  scope  limitations. 

The  remainder  of  this  chapter  highlights  the  literature  review  conducted  in  this 
field  as  well  as  the  goals  of  this  research.  Chapter  II  lays  out  the  simulation  methodology 
and  the  selected  guidance  laws  in  detail.  Chapter  III  describes  the  experimental 
procedures,  results  and  analysis,  while  Chapter  IV  presents  the  research  conclusions  and 
suggestions  for  further  research. 
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A. 


MISSILE  GUIDANCE  LAWS  LITERATURE  REVIEW 


Simply  stated,  the  goal  of  guidance  is  to  reach  the  target  [7].  In  general,  missile 
command  guidance  can  occur  in  two  basic  forms.  The  first  form  is  homing  guidance 
where  the  missile  relies  on  its  onboard  seeker  to  detect  the  target  and  compute  the 
required  command  guidance  through  integral  software  logic  circuits.  The  second  form 
relies  on  an  external  source  that  detects  the  target  and  the  missile,  computes  the  required 
guidance  to  bring  the  missile  towards  the  target  and  then  transmits  that  information  to  the 
missile  for  flight  control.  Such  command  guidance  control  is  advantageous,  as  the  missile 
is  not  required  to  have  a  seeker  onboard,  which  is  a  costly  component.  Such  guidance 
tends  to  result  in  very  good  missile  perfonnance.  However,  command  guidance  is  highly 
susceptible  to  tracking  errors.  The  quality  of  the  commanded  guidance  is  only  as  good  as 
the  quality  of  the  tracking  data.  Since  the  external  source  usually  remains  relatively 
stationary  in  space,  the  intercept  between  the  missile  and  the  target  usually  occurs  far 
away  from  the  command  source.  Hence  measurement  accuracy  of  the  tracking  data  and 
guidance  accuracy  are  limited. 

Conversely,  for  homing  guidance,  having  a  seeker  onboard  means  additional  cost, 
but  at  the  same  time  delivers  improved  guidance  accuracy  results  since  the  seeker  is 
continuously  approaching  the  target  as  time  progresses.  It  is  thus  apparent  that  some  fonn 
of  terminal  homing  guidance  is  highly  desirable  for  Air  to  Air  Missiles  (AAM),  which  are 
usually  engaging  targets  at  some  distance  from  the  launching  platfonn  and  are  highly 
maneuverable  at  the  same  time. 

In  general,  regardless  of  command  or  homing  guidance,  a  guidance  law  ultimately 
acts  as  the  detenninant  on  how  a  particular  set  of  commands  for  guidance  is  to  be 
generated.  In  general,  Goodstein  [8]  describes  three  guidance  law  general  categories  into 
which  all  other  guidance  laws  can  be  categorized.  There  are  many  special  cases  that  are 
modifications  of  these  three  basic  guidance  concepts  (see  Figure  1): 

•  LOS 

•  Pursuit 

•  Proportional 
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GUIDANCE  LAW  TYPES 


Figure  1.  Graphical  representation  of  three  basic  guidance  laws.  After  [8]. 

For  homing  guidance  implementations,  only  pursuit  guidance  and  proportional 
guidance  are  applicable.  LOS  guidance  by  definition  requires  an  external  control  vehicle 
to  establish  the  LOS  between  the  target  and  the  vehicle  along  which  the  missile  is  then 
guided  to  travel.  Hence,  for  this  study  into  advanced  guidance  laws,  the  focus  will  remain 
on  onboard  seeker-capable  implementations. 

1.  Pursuit  Guidance 

As  described  earlier,  pursuit  guidance  works  by  aiming  the  missile  directly  at  the 
target  throughout  the  entire  engagement.  It  is  a  simple  implementation  that  is  less 
sensitive  to  noise.  However,  it  is  not  effective  against  moving  targets  like  aircraft,  as  it 
usually  ends  up  in  an  energy-consuming  tail  chase  scenario.  There  are  other  applications 
where  the  speed  advantage  of  the  interceptor  is  very  large  (as  in  an  air-to-surface 
engagement  against  a  fixed  target),  in  which  case  pursuit  guidance  is  an  effective 
guidance  law. 
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2. 


Lead  Pursuit  /  Lead  Collision  Guidance 


Variations  of  the  pursuit  guidance  law  include  lead  pursuit  and  lead  collision  (see 
Figure  2).  As  the  name  implies,  lead  pursuit  means  that  the  missile  is  flying  a  course 
whereby  it  is  in  pursuit  of  a  leading  point  just  slightly  ahead  of  the  target.  As  this 
guidance  law  aims  to  predict  slightly  ahead  of  where  the  target’s  next  position  will  be,  it 
tends  to  be  more  effective  than  pure  pursuit  guidance  and  is  usually  able  to  engage  targets 
earlier  in  the  flight  path.  However,  it  still  has  essentially  the  same  problems  with  pure 
pursuit  guidance  and  is  seldom  used  in  systems  that  require  intercept  of  high  speed  high 
maneuver  targets. 

Lead  collision  is  a  further  extension  of  the  lead  pursuit  guidance.  It  is  also 
potentially  the  most  efficient  and  optimal  missile  trajectory  as  it  basically  involves 
pointing  the  missile  at  the  point  ahead  of  the  target  where  a  collision  would  occur  if  the 
target  continued  in  a  straight  line  with  no  acceleration.  This  is  an  efficient  trajectory,  as  it 
requires  minimal  control  effort  from  the  missile.  The  largest  drawback  is  that  it  requires 
the  target  to  fly  a  constant  trajectory  with  minimal  accelerations. 


TARGET 


INTERCEPT  POINTS 


PP 


Figure  2.  Pursuit  guidance  trajectories.  From  [9]. 
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3.  Proportional  Navigation 


As  one  of  the  most  popular  guidance  laws,  PN  or  some  augmented  form  exists  in 
many  missile  systems  in  the  world  today.  It  is  based  largely  upon  the  instantaneous 
direction  of  the  target  relative  to  the  missile  and  its  first  derivative  with  respect  to  time. 
There  are  two  generic  classes  of  PN,  pure  and  true.  Pure  PN  applies  the  commanded 
acceleration  with  reference  to  the  velocity  vector  of  the  missile;  whereas  True  PN  applies 
the  commanded  accelerations  with  reference  to  the  LOS  (see  Figure  3).  PN  has  a  highly 
nonlinear  set  of  governing  equations,  and  attempts  to  solve  them  tend  to  take  the 
approach  of  true  PN,  which  is  more  mathematically  tractable  as  compared  to  pure  PN. 
However,  in  practical  application  terms,  pure  PN  is  the  more  natural  PN  law,  as 
implementing  an  acceleration  vector  perpendicular  to  the  LOS  as  required  by  true  PN  is  a 
physical  challenge  in  practical  applications. 


PURE  PROPORTIONAL  NAVIGATION  TRUE  PROPORTIONAL  NAVIGATION 

Figure  3.  Geometry  showing  application  of  acceleration  vector.  From  [10]. 
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B.  GOALS 

This  research  thesis  is  focused  on  two  areas:  1)  to  investigate  the  perfonnance  of 
advanced  guidance  laws  (APN  and  DG)  and  compare  them  with  baseline  PN 
perfonnance,  2)  to  investigate  the  effects  of  noise  injection  into  the  seeker  and  analyze 
the  degradation  in  performance  between  the  different  compared  guidance  laws. 

This  will  be  achieved  by  modeling  the  system  utilizing  MATLAB®  Simulink® 
and  six  degrees  of  freedom  (6  DOF)  models.  The  efficacies  of  the  guidance  laws  will  be 
compared  using  the  kinematic  boundary  concept  espoused  by  Broadston  in  an  earlier 
paper  [11].  The  kinematic  boundary  is  basically  a  visual  representation  of  the  engagement 
envelope  with  the  target  at  the  center,  and  the  boundary  represents  the  maximum  range  at 
which  a  missile  can  be  launched  at  the  target  and  be  expected  to  hit.  The  missile  is 
assumed  to  be  launched  directly  at  the  target  at  all  heading  angles  and  it  undergoes  a 
short  period  of  constant  thrust  to  achieve  maximum  speed.  The  missile  subsequently 
glides  to  the  target  while  slowing  down  due  to  drag  forces.  The  target  is  assumed  to 
maintain  a  constant  speed. 

For  the  noise  study,  defined  baseline  noise  is  added  to  the  seeker  system  as  band 
limited  white  noise,  and  the  scenarios  are  run  over  a  hundred  simulations  to  determine  the 
percentage  of  hits  at  various  factors  of  the  baseline  noise. 
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II.  SIMULATION  METHODOLOGY 


For  the  simulation  model,  Broadston’s  work  on  a  simplified  6DOF  model  was 
utilized  as  the  main  simulation  engine  [11].  Part  of  this  thesis  investigates  the 
assumptions  that  were  utilized  in  creating  the  original  model  and  the  amendments  that 
were  made  to  update  the  model  for  use  in  this  thesis  paper. 

A.  SUMMARY  OF  SIMULATION  MODEL 

1.  SIX  DEGREES  OF  FREEDOM  (6DOF)  EQUATIONS 

For  a  free  body  in  space,  it  is  possible  to  be  translated  and  rotated  along  the  three 
principal  axes.  These  six  freedoms  of  movement  are  the  6DOF  that  an  unconstrained 
body  in  free  space  will  be  able  to  experience. 

The  general  practice  for  an  aircraft  body-centered  (ABC)  coordinate  frame  is  to 
define  the  origin  at  the  center  of  gravity  (c.g.)  of  the  body  with  the  x-axis  pointing 
towards  the  nose,  the  y-axis  pointing  towards  the  right  wing  (looking  at  the  aircraft  from 
top  down,  nose  pointing  up)  and  the  z-axis  is  90°  to  both  x  and  y  axes  pointing  straight 
down. 

For  simulation  of  the  air-to-air  combat  scenarios,  the  North-East-Down  (NED) 
frame  is  used  as  the  reference  frame.  The  NED  frame  has  its  origin  point  placed  at  the 
earth’s  surface,  with  the  x-axis  pointing  due  north,  y-axis  pointing  due  east  and  the  z-axis 
is  pointing  down  towards  the  center  of  the  earth. 

For  this  simulation,  flat  earth  approximations  are  used,  as  the  ranges  involved  in 
air-to-air  combat  are  relatively  short  as  compared  to  tactical  ballistic  missiles,  and  a  point 
mass  model  will  be  assumed  for  the  flight  dynamics.  The  following  vector  equations  fully 
describe  the  motion  dynamics  of  a  free  body  in  space  [12],  [13]: 
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(oB  =  -  J  lClBJcoB  +  J  lTB 
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Here  vB,(oB,g0,pNED,FB,TB, 


(force  -  translational  kinematics) 

(moment  -  rotational  dynamics)  (2  1) 

(attitude  -  rotational  kinematics) 

(navigation  -  translational  kinematics) 
are  vectors  and  Qb,C1,J,Bb,  are  square 


matrices.  The  equations  used  in  (2.1)  are  intended  for  flat  earth  approximations.  Hence 
they  do  not  include  terms  that  transfonn  the  NED  frame  to  an  earth  centered  inertial 
(ECI)  frame.  Those  tenns  would  be  needed  for  simulations  that  require  modeling  of  the 
missile  flight  path  over  a  spherical,  rotating  earth,  such  as  simulations  of  ballistic  missile 
trajectories. 


It  is  important  to  note  that  the  equations  in  (2.1)  provide  no  direct  correlations 
between  the  force  and  moment  equations.  In  the  physical  world,  a  dynamically  stable 
flight  body  would  have  its  c.g.  forward  of  its  center  of  pressure  (c.p.).  This  would  result 
in  the  missile  self-aligning  itself  to  the  relative  wind  direction  during  flight.  Therefore,  to 
simulate  this  stable  dynamic  behavior,  a  proportional-differential  controller  is  designed  to 
model  the  missile  attitude  such  that  it  simulates  actual  physical  behavior. 


2.  MISSILE  MODEL 


In  order  not  to  duplicate  previous  work,  this  thesis  adopts  Broadston’s  AIM- 120 
AMRAAM  model  [11]  with  its  flight  characteristics  and  dynamics.  The  missile  model 
was  created  based  on  capabilities  reported  in  open  source  literature  and  on  engineering 
approximations.  Hence  it  is  not  meant  to  be  an  exact  replica  of  the  actual  missile 
capabilities.  However,  the  simulation  model  is  created  modularly  so  the  missile  model 
characteristics  can  be  easily  modified  and  inserted  into  the  simulation  as  required.  The 
missile  body  dimensions  used  in  this  simulation  is  given  in  Table  1  and  has  been 
simplified  to  follow  the  models  described  in  Blakelock  [14]  and  Zarchan  [3]. 
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DESCRIPTION 

SYMBOL  USED  IN  CODE 

VALUE 

MISSILE  BODY  DIMENSIONS 

Missile  mass 

MASS 

156.8  kg 

Missile  diameter 

DIAM 

0.1778  m 

Missile  length 

LENGTH 

3.657  m 

Location  of  c.g.  measured  from 

the  nose 

XCG 

1.8288  m 

Length  of  the  nose  cone 

LN 

0.6769  m 

MISSILE  TAILPLANE  DIMENSIONS 

Hinge  line  distance  from  nose 

XHL 

3.454  m 

Tail  root  chord 

CRT 

0.4061  m 

Tail  tip  chord 

CTT 

0.0676  m 

Tail  extension 

TXT 

0.0676  m 

Tail  height 

HT 

0.2286  m 

MISSILE  WINGPLANE  DIMENSIONS 

Wing  to  radome  tangency  point 

distance  from  nose 

XW 

1.134  m 

Wing  root  chord 

CRW 

0.3554  m 

Wing  tip  chord 

CTW 

0  m 

Wing  extension 

WXT 

0  m 

Wing  height 

HW 

0.1778  m 

Table  1.  Summary  of  missile  dimensions  used  in  code.  From  [11]. 
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a.  Thrust  Characteristics 

The  AMRAAM  thrust  profile  is  simulated  as  a  constant  23,000  N  thrust 
for  the  initial  six  seconds  of  the  simulation  that  accelerates  the  missile  up  to  a  maximum 
speed  of  around  1,100  m/s,  which  is  approximately  Mach  3.5  at  an  altitude  of  6,000  m. 

b.  Moments  of  Inertia 

The  missile  is  modeled  as  a  thin  rod  for  the  y  and  z  axes  because  the 
missile  control  surfaces  are  assumed  to  have  minimal  impact  on  the  moment  of  inertias 
about  the  axes.  The  cylindrical  model  was  selected  for  the  x  axis. 

c.  Drag  Model 

In  general,  the  drag  force  is  a  combination  of  friction  drag  and  drag  caused 
when  the  integral  of  pressure  over  the  whole  surface  area  of  the  missile  body  is  nonzero 
[13].  The  components  of  total  drag  on  the  body  can  be  broken  down  into  three  major 
components,  parasitic  drag  (which  includes  friction  drag  and  form  drag),  induced  drag 
and  wave  drag.  In  reality,  the  three  drag  components  are  not  independent  and  cannot  be 
linearly  added  to  derive  total  drag.  However,  for  the  assumptions  of  the  simulation,  a 
simplified  drag  model  is  derived  from  the  first  two  components,  parasitic  and  induced 
drag.  The  following  equation  describes  the  total  drag  force  along  the  x-axis  of  the  ABC 
frame  [15]: 

f  =  (2.2) 

Parasitic  drag,  Cdo,  is  estimated  using  typical  values  from  [12].  Figure  4 
plots  the  value  of  Cdo  at  various  Mach  numbers.  The  plot  shows  the  difference  in  parasitic 
drag  values  caused  by  the  presence  or  absence  of  the  rocket  thrust  plume  during  the  boost 
phase  and  the  gliding  phase. 
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Parasitic  Drag  Coefficients  Mach  Number 


Figure  4.  Parasitic  drag  coefficient  variations  with  Mach  no. 


For  induced  drag,  Cdi,  it  is  normally  estimated  as  a  function  of  angle  of 
attack.  In  this  simulation,  the  normal  forces  acting  on  the  missile  body  are  used  instead. 
For  subsonic  flight,  a  crude  approximation  of  Cdi  is  made  as  the  applied  normal  force 
times  the  maximum  value  of  Cdo  in  subsonic  flight.  For  supersonic  flight,  Cdi  is  computed 
using  the  following  [15]: 


Ca  = 


c 


N 


ne(AR) 


Where  the  normal  force  coefficient  is  given  by: 


(2.3) 


C 


F, 


N 


N 


pv2sref 


(2.4) 


For  the  simulation,  the  normal  forces  and  induced  drag  acting  on  the  y- 
and  z-  axes  are  computed  separately  before  being  summed  together  as  a  vector  to  derive 


the  total  induced  drag  acting  upon  the  system. 
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d.  Drag  Model  Validation 


As  the  figures  used  in  generating  the  drag  forces  are  estimates,  a  sanity 
check  was  done  to  ensure  that  the  drag  profile  of  the  missile  is  not  too  far  from  reality.  To 
do  that,  open  source  data  for  AIM- 120  D  claiming  a  max  range  of  approximately  72  km 
was  compared  against  the  simulation  model  running  the  drag  profile.  In  a  tail  chase 
scenario  without  target  maneuvers,  the  missile  was  able  to  achieve  a  five  meter  radius 
impact  on  the  target  at  a  max  range  of  76.6  km.  This  represents  an  error  of  about  six 
percent,  which  is  acceptable  in  validating  the  parasitic  drag  coefficient  numbers  and 
induced  drag  model  used  in  this  simulation. 


3.  Target  Dynamics  Model 


For  the  simulation,  the  target  dynamics  were  modeled  as  a  point  mass  model  and 
implemented  with  the  following  vector  equation  [12]: 


[0  1  0  0  0  OYxl  [O' 

x  0  0  0  -co  0  0  x  0 

y  0  0  0  1  0  0  y  0 

y  0  0  0  0  0  j  0 

z  OOOOOlz  0 

z  OOOOOOz  a_ 


(2.5) 


In  this  target  dynamic  model,  the  target  is  modeled  with  lateral  accelerations  that 
are  given  as  turn  rate,  co.  The  vertical  acceleration  is  given  as  a  direct  input,  az,  to  the 
subsystem,  which  is  equal  to  the  gravitational  acceleration. 


4,  Guidance  Law  Implementation 

For  the  simulation,  the  guidance  laws  are  invoked  through  MATLAB®  function 
calls,  which  compute  the  required  acceleration  commands,  ny  and  nz  for  the  horizontal 
and  vertical  plane.  To  assist  the  functions  in  deriving  the  acceleration  commands,  the 
model  provides  the  function  with  data  from  the  seeker  head  and  the  inertial  measurement 
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unit  (IMU)  processor.  The  seeker  computes  target  range,  range  rate,  vertical  and 
horizontal  azimuth,  as  well  as  azimuth  angular  rates,  while  the  IMU  processes  the 
missile’s  own  state  vectors  and  orientation. 

Some  advance  guidance  laws  require  the  entire  target  state  vector  as  additional 
inputs,  and  a  simple  alpha-beta-gamma  (abg)  filter  was  implemented  to  estimate  the 
target  states  for  the  guidance  law  inputs. 

Section  C  of  this  chapter  will  discuss  the  three  guidance  laws  that  were  chosen  to 
be  implemented  in  greater  detail. 

5.  Navigation  Model 

The  IMU  provides  navigational  data  that  accelerometers  and  gyros  would  provide 
in  real  life.  The  missile  orientation  in  the  fonn  of  Euler  angles,  velocity,  acceleration, 
angles  of  attack  and  body  rotation  rates  are  computed  by  the  IMU.  The  IMU  is  assumed 
to  provide  this  information  as  truth. 

6.  Noise  Model 

For  the  noise  model,  it  was  necessary  to  change  the  simulation  solver  from  a 
variable  step  ordinary  differential  equation  (ode)  45  (Donnand-Prince)  solver  to  a  fixed 
step  size  ode4  (Runge-Kutta)  solver.  This  is  to  ensure  that  the  noise  input  process  is 
constant  at  each  time  step  of  the  simulation.  This  results  in  slightly  less  accurate 
simulation  runs  as  the  fixed  step  size  has  to  be  set  at  a  large  enough  value  that  ensures 
that  each  simulation  run  does  not  take  too  long  while  at  the  same  time  not  so  large  that 
the  simulation  results  become  very  inaccurate.  For  the  noise  simulations,  the  selected 
value  for  the  fixed  step  size  was  0.01s  and  it  resulted  in  a  difference  in  the  results  of 
around  1%  when  compared  to  the  ode45  solver. 

To  model  noise,  band-limited  white  noise  was  added  to  the  seeker  measurement 
outputs  of  range,  range  rate,  FOS,  and  FOS  rate.  Due  to  some  issues  with  the  abg-filter 
simulation  results,  it  was  decided  that  the  noise  for  the  filter  outputs  would  be  modeled  in 
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a  similar  fashion  to  the  seeker  output,  where  the  band-limited  white  noise  was  added 
directly  to  the  measurement  outputs  of  target  position,  velocity  and  acceleration. 

The  baseline  noise  model  is  defined  in  Table  2,  where  the  accuracy  of  the  radar 
systems  in  defining  range  was  estimated  at  ±10  m,  and  range  rate  at  ±1  m/s.  The  accuracy 
of  LOS  and  LOS  rate  was  defined  as  ±1  mrad  for  both  measurements.  To  test  for 
increased  noise  in  the  system,  a  common  gain  factor  is  applied  to  all  power  spectral 
density  values  as  defined  in  the  baseline  noise  model.  The  multiple  of  the  gain  factor  then 
determines  the  amount  of  noise  that  the  system  is  able  to  withstand. 


Noise 

Parameters 

Range 

(m) 

Range  Rate 
(m/s) 

LOS 

(rad) 

LOS  Rate 
(rad/s) 

Target 

Position 

(m) 

Target 

Velocity 

(m/s) 

Target 

Acceleration 

(m/s2) 

Power  Spectral 
Density 

1 

le-2 

le-8 

le-8 

1 

le-2 

2e-2 

Standard 

Deviation 

9.917 

0.992 

9.9e-4 

9.9e-4 

9.917 

0.992 

1.402 

Variance 

98.34 

0.983 

9.8e-7 

9.8e-7 

98.34 

0.983 

1.967 

Table  2.  Applied  baseline  noise  values  for  noise  simulation. 


B.  MODEL  ANALYSIS  AND  MODILICATIONS 

The  simplified  6DOF  model  that  was  adopted  for  use  in  this  thesis  was  modified 
from  the  original  in  order  to  solve  some  anomaly  issues  that  were  found  while  attempting 
to  generate  the  kinematic  boundary  comparisons.  In  the  initial  implementation  of  the 
model,  it  was  found  that  certain  anomalies  kept  appearing  in  the  kinematic  boundary 
plots,  an  example  of  which  can  be  seen  in  Figure  5. 
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Figure  5.  KB  plot  of  PN  vs.  non-maneuvering  target  showing  anomaly. 

There  was  a  predominance  of  anomalies  happening  at  the  broadsides  and  towards 
the  rear  quadrant.  In  the  case  of  Figure  5,  the  anomaly  occurs  at  105°  heading,  with  a 
sudden  decrease  in  the  maximum  range.  To  investigate  this  problem,  multiple  runs  were 
conducted  at  the  105°  heading  at  100  m  range  step  increments  and  the  minimum  miss 
distance  at  each  range  was  recorded  and  plotted  in  Figure  6. 


Min  Miss  Distance  vs  Range  Fired  (6)  105  deg  for 
PN  non-maneuvering  target 


Range  Missile  is  fired  at  target  (m) 


Figure  6.  Minimum  miss  distances  vs.  range. 
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As  is  evident  in  the  plot,  there  seems  to  be  a  specific  range  bin  whereby  the 
missile  miss  distances  at  end  game  increases  suddenly  before  going  back  to  the  baseline 
miss  distance.  The  sharp  increase  in  miss  distances  at  the  far  right  side  of  the  plot 
represents  the  maximum  range  of  the  missile. 

Upon  investigation,  it  was  discovered  that  the  region  where  the  anomalies  occur 
actually  corresponds  to  the  point  where  the  missile  is  decelerating  through  the  transonic 
region.  This  occurs  at  the  range  where  the  missile  is  fired  such  that  the  missile  is  at  the 
end  game  near  the  target  just  as  its  velocity  traverses  the  transonic  region.  Below  Mach 
one,  the  previous  drag  model  (see  Figure  7)  demonstrated  a  sudden  decrease  in  drag 
coefficient,  which  translates  to  lower  drag  forces.  The  guidance  law  perceives  this  as  a 
sudden  acceleration  on  the  part  of  the  missile,  and  simulation  studies  indicate  a  miss 
(where  a  miss  is  defined  as  a  minimum  miss  distance  greater  than  5  m).  It  can  be  seen 
that  if  the  missile  is  launched  even  further  away  such  that  it  passes  through  the  Mach  one 
boundary  before  entering  end  game,  the  missile  is  actually  able  to  impact  the  target  again. 

1.  Modified  Parasitic  Drag  Curve 

In  a  bid  to  reduce  the  transonic  region’s  instability,  the  drag  mode  as  proposed  by 
Broadston  [11]  was  modified  slightly  to  allow  for  a  gradual  increase  and  decrease  in  drag 
coefficient  through  the  transonic  boundary  instead  of  as  a  sudden  step  increase.  The  old 
and  new  parasitic  drag  models  are  shown  in  Figure  7  as  a  comparison. 

Original  Parasitic  Drag  Coefficient  vs  Mach  Number  Modified  Parasitic  Drag  Coefficient  vs  Mach  Number 


Figure  7.  Original  (left)  and  modified  (right)  induced  drag  models. 
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2. 


Additional  Fx  Saturation 


The  total  maximum  g-force  the  missile  is  assumed  to  be  able  to  sustain  is  52  times 
the  nonnal  gravitational  acceleration.  An  attempt  to  limit  the  commanded  acceleration  on 
all  three  axes  proportionally  led  to  algebraic  errors  in  the  simulation  runs.  Therefore,  a 
simplified  method  was  adopted  to  limit  the  accelerations,  which  basically  is  a  direct 
limitation  placed  on  the  accelerations  in  the  individual  axes  not  to  exceed  30g  each. 
Initially,  only  the  commanded  accelerations  in  the  y-  and  z-  axes  were  limited.  However, 
in  the  simulation  analysis  and  verification  study,  it  was  discovered  that  the  x-  axis 
acceleration  must  be  limited  as  well  to  provide  consistent  and  smooth  results.  Therefore 
the  final  model  contains  limitations  on  all  three  axes.  The  limiters  for  the  x-,  y-  and  z- 
axis  are  modeled  as  saturation  blocks  within  the  Simulink®  model. 

3.  Additional  Turn  Limiter 

In  the  simulation  analysis  stage  of  the  project,  it  was  discovered  that  using  the 
Simulink®  switch  block  alone  to  govern  the  moment  when  the  target  is  expected  to 
execute  a  constant  g-turn  maneuver  towards  the  missile  could  lead  to  unexpected  effects 
in  certain  simulation  conditions.  The  switch  block  is  basically  given  the  time-to-go  value 
(as  derived  from  range-to-go  and  range -to-go  rate)  and  once  the  time-to-go  drops  below 
the  simulation  detennined  constant  value,  it  passes  the  constant  turn  command  through  to 
the  target  dynamics  model.  It  was  discovered  that  in  certain  end  game  simulations  where 
the  missile  is  in  the  transonic  region  of  its  flight,  the  sudden  decrease  in  the  missile’s 
velocity  caused  the  time-to-go  value  to  increase  back  above  the  determined  constant 
value.  This  causes  the  switch  block  to  stop  passing  the  turn  command.  To  prevent  such 
situations  from  occurring  and  introducing  added  anomalies  to  the  results,  a  turn  limiter 
was  coded  into  the  simulation.  The  turn  limiter  basically  holds  the  turn  command  once  it 
is  passed  through  the  switch  block. 
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c. 


GUIDANCE  LAWS 


This  section  describes  the  three  guidance  laws  that  were  selected  for  investigation 
in  detail.  The  first  two  laws,  PN  and  augmented  proportional  navigation  (APN)  were 
selected  as  the  two  best  performing  guidance  law  that  Broadston  [11]  investigated.  From 
the  additional  literature  research,  an  additional  guidance  law  that  was  not  investigated 
previously  using  the  kinematic  boundary  technique  was  selected  to  compare  against  the 
other  two  laws.  This  additional  law  is  based  on  differential  geometry  (DG)  [4],  [5],  [16], 
and  [17],  and  basically  generates  commanded  accelerations  as  a  function  of  the  missile 
flight  geometry.  The  geometry  of  a  typical  engagement  scenario  is  shown  in  Figure  8. 


Figure  8.  Missile  engagement  geometry.  After  [11]. 


1.  Proportional  Navigation  (PN) 

As  described  briefly  in  the  literature  review  section  of  the  introduction  chapter, 
the  implementation  of  the  PN  guidance  law  in  this  simulation  is  modeled  after  the  True 
PN  model  where  the  commanded  accelerations  are  applied  perpendicularly  to  the  LOS. 
The  2D  commanded  acceleration  is  given  as  follows  [3]: 

nc  =  N'VC&  (2.6) 
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While  it  is  mathematically  more  tractable  to  compute  the  commanded  acceleration 
perpendicular  to  the  LOS,  in  the  simulation,  it  is  instead  easier  to  apply  the  commanded 
acceleration  perpendicular  to  the  missile  body  axis  instead.  In  order  to  achieve  this,  the 
commanded  acceleration  obtained  in  (2.6)  is  corrected  to  the  missile  y-  axis  as  follows: 


,  N’VC& 
cosier,) 


(2.7) 


The  corrected  missile  acceleration  is  then  fed  directly  to  the  missile  dynamics 
block  for  simulation.  The  navigational  constant,  N’  is  selected  to  be  five  for  the  entire 
simulation  for  PN  law.  As  described  in  Zarchan  [3],  setting  the  navigation  ratio  at  five 
reduces  the  required  missile  acceleration  advantage  to  1.67,  which  increases  the 
maximum  distance  that  the  guidance  law  will  be  able  to  hit  the  target. 


2.  Augmented  Proportional  Navigation  (APN) 


An  advanced  version  of  PN  law,  the  APN  law  is  basically  PN  augmented  with  a 
component  that  accounts  for  the  maneuvering  target  dynamics.  As  described  in  [3],  an 
alternative,  mathematically  equivalent,  expression  for  (2.6)  is  as  follows: 


Where  tg0  is  the  time  to  go  until  intercept  and  yr  is  defined  in  Figure  8.  The  terms 
in  the  parentheses  are  also  referred  to  as  the  zero  effort  miss  (ZEM)  [3].  ZEM  basically 
states  the  miss  distance  that  would  result  if  the  missile  and  the  target  experienced  no 
further  accelerations  and  were  allowed  to  keep  traveling  with  the  same  velocity. 
Therefore,  if  we  were  to  assume  a  target  that  maneuvers,  than  the  ZEM  equation  must  be 
augmented  with  an  additional  term  that  addresses  the  maneuver.  The  new  equation  for 
ZEM  with  target  maneuver  is  as  follows: 


ZEM 


new 


v  + v t  h — nr 

-  r  Sr  go  2  ‘  S 


(2.9) 


Where  nt  is  the  target  acceleration.  Therefore,  the  augmented  PN  law  becomes  as 
follows: 
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(2.10) 


. .  N n. 
n„  =  N  Vr.cr  + - -  = 


N'(yr  +  yrtgo  +  -n,t2g0 ) 


t 


go 


It  can  thus  be  seen  that  APN  is  basically  PN  law  augmented  with  a  constant 
proportion  of  the  target  acceleration.  This  increases  the  complexity  of  the  system  as  a 
tracking  filter  that  will  estimate  the  target  acceleration  is  now  required. 


The  implementation  of  APN  in  the  simulation  is  given  by  the  vector  equation 
from  [18]: 
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(2.11) 


According  to  [3],  the  optimal  value  for  N’  occurs  when  N’  =  3.  This  would  be 
evaluated  subsequently  in  the  scenario  analysis.  From  this  vector,  only  the  values  for  y- 
and  z-  accelerations  are  used  as  x-  is  computed  from  thrust  and  drag  components  of  the 
model. 


3.  Differential  Geometry  (DG) 

DG  guidance  is  a  form  of  curvature  control  command  where  the  command  for 
missile  guidance  is  developed  based  on  the  curvature  concept  in  curve  theory  and  the 
relative  rotational  motion  of  a  pseudo-missile  pointing  velocity  vector  [16]. 

According  to  Chiou  in  [4],  the  Frenet-Serret  formula  for  classical  differential 
geometry  curve  theory  was  utilized  in  developing  this  control  law.  The  fundamental 
theory  lying  behind  PN  and  APN  guidance  laws  is  that  a  fixed  model  for  the  likely  target 
maneuvers  have  to  be  assumed  and  is  then  used  to  derive  the  guidance  law.  This 
approach  has  certain  drawbacks,  not  least  of  which  is  that  the  real  world  target  seldom 
exhibit  constant  maneuver  tactics.  Therefore,  the  DG  guidance  law  attempts  to  address 
this  issue  by  ignoring  specific  target  maneuver  models  in  its  formulation. 
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Figure  9.  Geometry  of  engagement  scenario.  From  [17]. 


A  simple  2D  application  of  the  DG  guidance  law  in  the  time  domain  is  described 
in  [17],  where  the  engagement  geometry  is  laid  out  as  in  Figure  9. 

rq 


K.=N2Kt^^-  +  N'  2 


(2.12) 


c°s r/m  K,  cos//,,, 

In  this  fonnulation,  Kt  and  Km  represents  a  measure  of  the  curvature  of  the  target 
and  the  missile’s  trajectories  respectively.  r|t  and  r|m  represents  the  lead  angles  of  the 
target  and  the  missile,  r’  represents  the  closing  speed  of  the  system  along  the  LOS  and  q’ 
represents  the  angular  rate  of  LOS.  N’  is  the  control  gain  and  is  not  to  be  confused  with 
N,  which  is  defined  as  the  ratio  of  the  magnitude  of  the  target  and  missile  velocities,  vt 
and  vm. 


N  = 


The  curvature  of  a  circular  track  can  be  described  simply  as: 

Ml 


Substituting  equation  (2.13)  and  (2.14)  into  (2.12): 
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23 


(2.16) 


cos  TJ.  .  Vo 

- —  +  N  — 2 — 

COS  77  cos  77, „ 

im  im 

From  the  respective  definitions,  cos(r|m)  =  cos(ol),  the  DG  law  expressed  in  (2.16) 
can  be  seen  to  be  the  basic  PN  law  augmented  by  a  geometric  tenn  in  front  that  requires 
target  acceleration  as  well  as  target  velocity  heading  information. 

For  implementation  within  the  simulation,  the  acceleration  requirements  are 
computed  respectively  in  the  horizontal  and  the  vertical  plane  to  generate  the  required  y- 
and  z-  axes  acceleration. 
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III.  SIMULATION  SCENARIOS,  COMPARISON  AND  ANALYSIS 


To  maintain  similarity  with  Broadston’s  [11]  work,  the  simulation  scenario 
parameters  were  kept  largely  the  same  so  that  effective  comparisons  could  be  made.  As 
before,  the  simulations  are  run  at  a  fixed  6000  m  altitude  with  the  target  and  missile 
initial  speed  at  Mach  0.83.  The  missile  is  always  considered  to  be  fired  pointing  directly 
at  the  target  at  all  heading  angles,  and  for  the  maneuvering  scenarios,  the  target  is  always 
assumed  to  pull  a  6g  turn  maneuver  towards  the  missile  at  tg0  =  3  s.  Table  3  summarizes 
the  various  test  scenarios. 


Test  No. 

Title 

Description 

Simulation  Settings 

1A 

APN  Navigational 
Constant 

To  test  for  N'  =  3  and  N'  =  5  for 

APN  law  and  comparethe  results. 

ode4S  (Dormand-Prince)  simulation.  6g  maneuvering 
target3  s  priorto  impact. 

IB 

DS  Navigational 

Constant 

To  test  for  N'  =  3,  N'  =  5  and  N'  =  7 
for  DG  law  and  comparethe 
results 

ode45  (Dormand-Prince)  simulation.  6g  maneuvering 
target 3  s  priorto  impact. 

2A 

Non-maneuvering 

target 

Testforthe  effectiveness  of  each 
guidance  law  against  constant 
velocity  target. 

ode45  (Dormand-Prince)  simulation.  No  target 
maneuvers.  N'  =  5  for  all  guidance  laws. 

2B 

Maneuvering  target 

Testforthe  effectiveness  of  each 
guidance  law  against  6g 
maneuvering  target. 

ode45  (Dormand-Prince)  simulation.  6g  maneuvering 
target.  N'  =  5  for  all  guidance  laws. 

3A 

Non-maneuvering 
targetwith  noise 

Compare  and  contrast 
effectiveness  of  each  guidance  law 
upon  addition  of  noise  against 
constantvelocity  target 

ode4  (Runge-Kutta)  simulation.  No  target  maneuvers. 

N'  =  5  for  all  laws.  Noise  factor  lx  baseline  noise. 

3B 

maneuvering  target 
with  noise 

Compare  and  contrast 
effectiveness  of  each  guidance  law 
upon  addition  of  noise  against  6g 
maneuvering  target 

ode4  (Runge-Kutta)  simulation.  6g  maneuvering 
target.  N'  =  5  for  all  laws.  Noise  factor  lx  baseline 
noise. 

4 

Noise  tolerance 
study 

Compare  noise  tolerance  levels  of 
each  guidance  law  at  45°  and  135° 
heading  with  100  simulation  runs 
each.  For  both  maneuvering  and 
non-maneuvering  scenarios. 

ode4  (Runge-Kutta)  simulation.  N'  =  5  for  all  laws.  100 
simulation  runs.  Noise  factor  increased  until  just 
before  more  than  30%  misses.  Tested  at  90%  of  max 
range  at  each  heading. 

Table  3.  Summary  of  simulation  test  scenarios. 


For  the  kinematic  boundary  simulations,  the  scenarios  were  run  with  an  initial 
10000  m  launch  range  for  the  missile  and  target  at  a  tail  chase  scenario  of  0°  heading. 
The  simulation  then  slowly  increments  the  launch  range  with  1000  m  steps,  then  100  m 
steps,  up  to  a  resolution  of  10  m  steps  to  determine  the  max  range  whereby  the  missile 
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first  misses  the  target  by  more  than  5  m  at  the  end  game.  Once  the  maximum  range  that 
the  missile  can  hit  the  target  is  detennined  at  that  bearing,  the  heading  is  then 
incremented  by  5°  before  the  next  simulation  run  to  determine  the  maximum  range  at  that 
bearing. 

For  each  simulation  run,  there  are  three  stopping  criteria  which  would  stop  the 
simulation  once  any  one  of  the  three  is  met.  The  first  stopping  condition  is  a  velocity  stop 
and  it  stops  the  simulation  once  the  missile  speed  falls  below  the  target  speed.  The  second 
condition  is  a  range  rate  stop  and  it  stops  the  simulation  when  the  rate  of  change  of  the 
distance  between  the  missile  and  the  target  becomes  positive.  The  last  condition  is  a 
range  stop  and  it  is  triggered  when  the  range  between  the  missile  and  the  target  is  less 
than  0.0001  m. 

The  simulation  is  run  from  0°  to  180°  bearing  and  the  maximum  launch  range 
results  are  all  stored  and  finally  plotted  in  a  polar  plot  as  the  kinematic  boundary  of  the 
simulation.  The  time  taken  for  one  such  simulation  run  varies  from  as  little  as  2.5  hours 
for  the  simple  PN  guidance  law  with  ode45  solver  up  to  8.5  hours  for  the  more 
complicated  DG  law  in  a  6g  target  maneuvering  scenario  running  the  same  variable  step 
solver. 

For  the  noise  scenarios,  using  the  ode4  fixed  step  solver,  the  simulation  times 
range  from  7.5  hours  to  10  hours  for  a  single  kinematic  boundary  simulation. 

When  using  the  ode45  variable  step  solver,  PN  law  running  with  a  simulation 
tolerance  of  le-3  and  minimum  step  size  of  le-3  and  maximum  step  size  at  auto  provides 
reasonably  accurate  results  in  a  short  time.  However,  for  the  more  complicated  APN  law, 
the  tolerance  and  minimum  step  size  of  the  ode45  solver  has  to  be  tightened  to  le-5  and 
the  maximum  step  size  specified  at  0.1  for  reasonably  accurate  results. 

For  the  DG  law,  even  tighter  simulation  parameters  are  required.  While  the 
tolerance  and  minimum  step  size  is  held  similar  to  APN,  the  maximum  step  size  has  to  be 
decreased  to  0.02  and  the  filter  sample  rate  for  the  abg-filter  set  at  50  Hz  for  reasonable 
results  to  be  generated.  In  contrast,  the  APN  law  was  able  to  generate  reasonable  results 
with  the  abg-filter  set  to  a  10  Hz  sample  rate  with  the  maximum  step  size  at  0.1. 


26 


For  the  noise  scenarios,  the  fixed  step  size  is  set  at  0.01  as  it  was  empirically 
tested  to  provide  the  most  accurate  results  for  all  3  guidance  laws  within  the  least  amount 
of  simulation  run  time. 

The  simulation  SIMULINK  models  are  shown  in  Appendix  A  while  the 
MATLAB®  codes  are  organized  within  Appendix  B. 

To  demonstrate  the  validity  of  the  simulation  model,  the  engagement  geometry 
and  time  history  plots  of  certain  parameters  are  shown  in  Appendix  C. 
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A. 


TEST  1A  -  APN  NAVIGATIONAL  CONSTANT 


In  this  test,  the  APN  law  was  tested  for  the  effects  of  the  values  chosen  for  the 
navigational  constant.  In  [8],  N’  =  3  was  derived  mathematically  as  the  optimal  value  for 
the  guidance  law.  The  results  for  this  scenario  are  shown  in  Figure  10.  The  polar  plot 
represents  the  kinematic  boundary  of  the  guidance  law  at  the  two  navigational  constant 
while  the  top  right  hand  plot  is  the  same  plot  plotted  linearly  to  emphasize  the  difference. 
The  bottom  right  hand  plot  directly  plots  the  difference  in  max  ranges  between  APN 
N’=3  and  APNN’=5. 


Bearing  of  launch  (deg)  0  cfeg  represents  tailchase 


Figure  10.  Comparison  of  APN  with  N’=3  and  N’=5. 

It  can  be  seen  from  the  difference  plot  that  there  is  a  significant  improvement  to 
the  performance  of  the  missile  when  N’=5  especially  between  40°  and  80°  bearing  where 
improvement  of  max  range  between  2000  m  and  8000  m  are  seen.  For  the  rest  of  the 
bearings,  there  is  a  minor  advantage  to  N’=3  but  the  difference  is  insignificant. 

It  is  clear  that  having  the  navigational  constant  set  at  N’=5  offers  a  general 
improvement  to  the  performance  of  the  missile  hence  for  APN  law,  N’  is  fixed  at  N’=5 
for  the  rest  of  the  simulations. 
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B. 


TEST  IB  -  DG  NAVIGATIONAL  CONSTANT 


As  in  test  1A,  the  DG  law  is  tested  at  various  navigational  constant  values  for 
maneuvering  target  scenario.  N’  =  3,  5  and  7  were  chosen  for  the  test  and  the  results  are 


shown  in  Figure  1 1 . 


Figure  11.  Comparison  of  DG  with  N’  =  3,  5  and  7. 


It  can  be  clearly  seen  that  N’=7  resulted  in  the  best  performance  envelope  of  the 
missile  in  this  maneuvering  scenario.  It  showed  significant  improvements  up  to  10000  m 
over  N’=3.  However,  it  demonstrates  only  a  slight  improvement  of  up  to  2000  m  over 
N’=5  and  only  in  a  very  narrow  tail  chase  region  of  approximately  0°  to  20°. 

Since  there  is  only  a  slight  difference  between  N’=5  and  N’=7,  the  navigational 
constant  for  DG  is  selected  to  be  N’=5  so  as  to  be  consistent  with  the  other  two  guidance 
laws  that  are  being  tested. 
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c. 


TEST  2A  -  NON-MANEUVERING  TARGET 


In  this  test,  each  of  the  three  guidance  laws  was  simulated  against  a  constant  0.83 
Mach  velocity  target.  As  both  APN  and  DG  laws  are  basically  PN  laws  augmented  with 
an  additional  control  term  that  takes  into  account  the  target  acceleration,  it  is  predicted 
that  all  three  guidance  laws  should  perform  similarly  when  the  target  exhibits  constant 
velocity  characteristics  for  the  entire  simulation.  Figure  12  summarizes  the  results  for  all 
three  guidance  laws. 


Non- maneuvering  comparison  for  PN,  APN  &  DG 
90  150000 


x  to1  Plot  of  Max  Ranges  for  PN,  APN  &  DG  (non-maneuver) 


Plot  of  Difference  between  Max  Ranges  for  PN  against  APN  &  DG 


Figure  12.  Non-maneuvering  comparison  of  PN,  APN  and  DG 


It  can  be  seen  that  the  simulation  results  are  similar  to  the  predicted  results.  The 
kinematic  boundary  plot  shows  that  all  three  laws  are  within  1%  -  3%  error  ranges  of 
each  other  and  that  the  difference  plot  shows  a  difference  of  -1000  m  to  3000  m  for  DG 
when  compared  to  PN  and  a  difference  of  up  to  1000m  for  a  narrow  heading  around  120° 
for  APN  when  compared  to  PN. 

This  difference  amounts  are  of  little  significance  and  it  can  be  concluded  that  all 
three  guidance  laws  perform  similarly  in  a  constant  velocity  target  scenario. 
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D. 


TEST  2B  -  MANEUVERING  TARGET 


Now  that  the  non-maneuvering  results  show  that  all  three  guidance  laws  are 
similar  at  the  baseline  non-maneuvering  scenario,  the  target  is  modeled  with  constant  6g 
acceleration  towards  the  missile  when  tgo  -  3  s.  Figure  13  summarizes  the  results. 


Maneuvering  comparison  for  PN,  APN  &  DG 

90  150000 


Plot  of  Max  Ranges  for  PN,  APN  &  DG  (6g  maneuver) 
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Plot  of  Difference  between  Max  Ranges  for  DG  against  PN  and  APN 


J-  q.  -l  -l  u  .i  i _ i _ 

0  20  40  00  80  100  120  140  100  180 

Bearing  of  launch  (d eg)  0  deq  repieserts  talchase 


Figure  13.  Maneuvering  comparison  of  PN,  APN  and  DG. 


It  can  be  seen  from  the  results  that  all  three  guidance  laws  perform  similarly  in  the 
forward  and  broadside  quadrants  of  the  kinematic  boundary  where  the  engagement 
geometry  varies  from  90°  to  head  on  with  the  target  at  180°. 

The  rear  quadrant,  however,  shows  significant  differences  between  the  three 
guidance  laws.  It  can  be  seen  that  PN  is  the  worst  performing  guidance  law  in  a  tail  chase 
maneuvering  target  scenario  except  for  a  narrow  bandwidth  around  60°  heading  where  it 
performs  slightly  better  than  APN.  APN  in  contrast  behaves  generally  better  than  PN 
between  20°  to  50°  heading  and,  other  than  the  slight  inferior  performance  compared  to 
PN,  APN  is  essentially  performing  similarly  to  PN  guidance. 

The  law  that  clearly  performs  better  is  DG  where  improvements  of  up  to  12000  m 
in  maximum  range  can  be  achieved  in  the  rear  quadrant.  However,  careful  analysis  shows 
that  DG  has  a  slight  performance  issue  in  the  front  quadrant  when  compared  to  the  other 
two  guidance  laws  even  though  the  difference  is  insignificant. 
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E. 


TEST  3A  -  NON-MANUEVERING  TARGET  WITH  NOISE 


For  this  test,  baseline  noise  was  added  to  the  system  and  a  single  simulation  to 
derive  the  kinematic  boundary  was  carried  out.  While  the  noise  process  is  random  and  the 
kinematic  boundary  simulation  will  output  a  different  result  each  time  it  is  ran,  the 
variations  from  one  degree  step  to  the  next  should  even  out  across  the  entire  180° 
boundary.  This  hypothesis  was  put  to  the  test  and  the  simulation  results  for  each  of  the 
three  guidance  laws  in  a  non-maneuvering  scenario  was  compared  with  their  own  results 
in  a  noiseless  simulation  and  summarized  in  Figures  14,  15  and  16. 


Noise  comparison  for  PN  (no  maneuver) 

90 

150000 


4 

x  10  Plot  of  Max  Ranges  for  PN  0  turn  with  &  w/o  noise 


Plot  of  Difference  between  Max  Ranges  for  PN  Otum  with  &  w/o  noise 


Figure  14.  Noise  comparison  for  PN  (no  maneuver) 


It  can  be  seen  that  the  PN  law  under  the  influence  of  baseline  noise  is  hardly 
affected  by  it  and  is  able  to  perfonn  just  as  well  as  the  scenario  where  noise  was  not 
added. 
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Noise  comparison  for  APN  (no  maneuver) 
90 

150000 


4 

x  io  Plot  of  Max  Ranges  for  APN  0  Urn  with  &  w/o  noise 


Plot  of  Difference  between  Max  Ranges  for  APN  0  turn  with  &  w/o  noise 


Figure  15.  Noise  comparison  for  APN  (no  maneuver) 


For  the  APN  law,  it  can  be  seen  that  it  is  more  affected  by  baseline  noise  in  the 
system  as  compared  to  PN.  There  is  an  increasing  difference  in  maximum  ranges, 
especially  towards  the  front  quadrant  where  the  difference  is  as  much  as  3000  m. 


Noise  comparison  for  DG  (no  maneuver) 

90 

150000 


x  iq4  Plot  of  Max  Ranges  for  DG  0  turn  with  &  w/o  noise 


Plot  of  Difference  between  Max  Ranges  for  DG  0  turn  with  &  w/o  noise 


Figure  16.  Noise  comparison  for  DG  (no  maneuver) 
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It  can  be  seen  that  DG  suffers  approximately  the  same  level  of  degradation  as 


compared  to  APN  under  the  influence  of  baseline  noise. 


Noise  comparison  for  PN,  APN  &  DG  (no  maneuver) 


150000 


Plot  of  Difference  between  Max  Ranges  for  6g  turn  wth  &  w/o  noise 


Figure  17.  Noise  comparison  for  PN,  APN  and  DG  (no  maneuver) 


When  the  three  guidance  laws  are  compared  against  each  other  under  the 
influence  of  noise,  it  can  be  seen  from  Figure  17  that  the  difference  between  the  two 
advanced  guidance  laws  and  PN  has  increased  by  about  2000  m  when  compared  with  the 
results  in  test  2A. 

This  is  consistent  with  the  expectation  that  the  advanced  guidance  laws  should 
suffer  greater  perfonnance  penalties  due  to  the  additional  complexity  of  the  guidance  law 
and  the  added  command  acceleration  terms. 
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F. 


TEST  3B  -  MANEUVERING  TARGET  WITH  NOISE 


Similarly  to  test  3A,  the  simulation  is  now  carried  out  with  the  target  maneuver 


scenario.  The  results  for  the  individual  guidance  laws  are  shown  in  Figures  18,  19  and  20. 


Noise  comparison  forPN  (6g  maneuver) 
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x  io  Plot  °f  Max  Ranges  for  PN  6  turn  with  &  w/o  noise 


Plot  of  Difference  between  Max  Ranges  for  PN  6  turn  with  &  w/o  noise 
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Figure  18.  Noise  comparison  for  PN  (6g  maneuver) 


Under  the  no  target  maneuver  noise  scenario,  PN  was  the  only  guidance  law  out 
of  the  3  that  did  not  suffer  perfonnance  degradation  under  the  effects  of  noise.  Now  that 
the  target  is  maneuvering,  PN  is  seen  to  suffer  the  same  amount  of  performance  penalty 
as  APN  and  DG  did  under  the  no  maneuver  scenario.  The  spike  at  the  55°  heading  point 
that  seems  to  show  a  region  where  the  PN  law  performed  better  with  noise  then  without 
noise  is  most  probably  an  anomaly. 
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Noise  comparison  for  APN  (6g  maneuver) 

go 
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x  IQ4  Plot  Max  Ranges  for  APN  6  tim  with  &  w/o  noise 


Plot  of  Difference  between  Max  Ranges  for  APN  6  turn  with  &  w/o  noise 


Figure  19.  Noise  comparison  for  APN  (6g  maneuver) 


As  with  the  results  for  PN,  it  can  be  seen  that  APN  suffers  much  greater 
perfonnance  degradation  around  the  60°  to  70°  heading  region.  While  PN  ‘suffered’  a 
perfonnance  increase,  APN  suffers  a  performance  decrease.  The  reason  for  this  increased 
perfonnance  anomaly  region  is  most  likely  attributable  to  the  transonic  region  effect. 


270  Bearing  of  launch  (deg)  0  deg  represents  tailchase 


Figure  20.  Noise  comparison  for  DG  (6g  maneuver) 
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DG  does  not  seem  to  suffer  as  much  performance  degradation  due  to  the  effects  of 
noise  as  much  as  the  other  two  guidance  laws.  It  can  be  seen  that  DG  suffers 
approximately  the  same  amount  of  perfonnance  loss  due  to  noise  regardless  of  the  target 
making  or  not  making  a  maneuver  in  the  scenario. 


Noise  comparison  for  PN,  APN  &  DG  (6g  maneuver) 
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Figure  2 1 .  Noise  comparison  for  PN,  APN  and  DG  (6g  maneuver) 


From  Figure  21,  it  can  be  seen  that  the  addition  of  noise  in  the  simulation  does  not 
affect  the  relative  performance  of  the  DG  law  when  compared  to  the  other  two  laws.  DG 
is  still  the  best  perfonning  guidance  law  under  target  maneuvering  scenarios. 
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G. 


TEST  4  -  NOISE  TOLERANCE  STUDY 


In  the  noise  tolerance  study  section,  each  simulation  data  point  was  run  for  a  total 
of  100  times  with  random  white  noise  throughout.  Due  to  the  time  taken  to  run  each  data 
point  100  times  with  a  fixed  step  ode4  solver  at  0.01  s  step  size  and  100  Hz  filter  sample 
rate,  the  simulations  were  run  only  at  the  45°  and  135°  heading. 

For  each  of  the  three  guidance  laws,  and  for  both  maneuvering  and  non¬ 
maneuvering  target  scenarios,  the  baseline  noiseless  maximum  range  at  the  two  bearings 
was  first  determined.  A  10%  reduction  of  the  maximum  range  was  taken  and  the 
simulations  were  then  run  at  that  range.  The  baseline  noise  in  the  system  is  then  increased 
by  a  constant  factor  and  the  factor  was  increased  until  no  less  than  70%  of  the  100  runs 
were  registered  as  successful  hits  on  target  (i.e.  for  the  missile  to  come  within  5  m  radius 
of  the  target). 

The  results  are  summarized  in  Table  4. 


Target 

Heading 

(deg) 

Target 

Maneuver  (g) 

Guidance 

Law 

Original 
Max  Range 
(m) 

10% 

Reduced 
Range  (m) 

Factor  of 
Baseline 

Noise 

Percentage 
Hit  (%) 

PN 

46,800 

42,120 

17. 7x 

81 

17. 8x 

65 

0 

APN 

46,790 

42,110 

7.2x 

76 

7.3x 

60 

DG 

47,220 

42,500 

16. 2x 

82 

45 

16. 3x 

65 

PN 

39,190 

35,270 

16. Ox 

75 

16. lx 

69 

6 

APN 

46,680 

42,010 

7.3x 

80 

7.4x 

54 

DG 

47,030 

42,320 

16. 8x 

75 

16. 9x 

55 

PN 

92,530 

83,270 

4.9x 

97 

5. Ox 

69 

0 

APN 

92,320 

83,090 

5.3x 

78 

5.4x 

49 

DG 

90,140 

81,120 

4.5x 

82 

135 

4.6x 

38 

PN 

88,490 

79,640 

6. lx 

76 

6.2x 

38 

6 

APN 

91,200 

82,080 

5.5x 

94 

5.6x 

66 

DG 

88,780 

79,900 

4.7x 

83 

4.8x 

29 

Table  4.  Summary  table  for  noise  tolerance  of  guidance  laws. 
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From  the  results,  it  can  be  seen  that  at  the  45°  heading,  both  PN  and  DG  laws 
exhibit  robust  tolerance  to  noise  and  can  take  noise  factors  around  16x  to  17x  of  baseline 
for  both  the  maneuvering  and  non-maneuvering  scenarios.  However,  it  is  significant  to 
note  that  the  maximum  ranges  for  PN  law  is  nearly  7  km  shorter  than  DG  and  APN  in  the 
maneuvering  scenario.  In  contrast,  APN  suffers  larger  performance  degradation  under  the 
effects  of  noise  and  has  a  noise  tolerance  of  only  around  7x  though  at  a  maximum  range 
similar  to  that  of  DG. 

At  the  135°  heading,  all  three  guidance  laws  exhibit  similar  noise  tolerance 
characteristics  of  about  5x  the  baseline  noise. 

As  can  be  seen  in  test  3B,  the  three  guidance  laws  suffered  similar  levels  of 
perfonnance  degradation  at  the  forward  quadrant  of  the  engagement  scenario.  For  the  rear 
quadrant,  it  could  be  seen  that  PN  and  DG  suffered  similar  perfonnance  degradation 
when  compared  to  their  respective  noiseless  scenario  perfonnances  but  APN  suffered 
additional  perfonnance  losses  at  the  rear  quadrant.  That  result  is  reflected  in  the  noise 
tolerance  study  at  the  45°  heading  where  it  can  be  seen  that  the  greater  the  perfonnance 
degradation  the  guidance  law  suffered,  the  lower  the  noise  tolerance  of  that  guidance  law 
at  that  particular  heading. 

The  full  data  and  histogram  plots  for  each  of  the  scenarios  described  in  Table  4, 
with  a  hit  percentage  greater  than  70%,  is  presented  in  Appendix  D. 
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IV  CONCLUSION 


A.  CONCLUSIONS 

The  developed  simplified  6DOF  simulation  model  proved  to  be  a  useful  tool  in 
the  modeling,  simulation  and  evaluation  of  the  different  guidance  laws.  Even  with  a  set  of 
simplifying  assumptions  about  physical  system  perfonnance  characteristics,  the 
simulation  results  matched  expectations. 

As  a  result  of  the  simulations,  it  was  proved  that  a  navigational  constant  of  N’  =  5 
is  optimal  for  both  the  APN  law  as  well  as  the  DG  law. 

For  the  maneuvering  target  scenario,  this  study  has  shown  that  while  the  APN  law 
does  indeed  perfonn  better  than  the  PN  law  except  for  a  very  narrow  region,  the  new  DG 
guidance  law,  which  is  based  upon  target  motion  geometries  instead  of  expected  target 
maneuvering  models,  performed  the  best  out  of  the  three  guidance  laws  tested.  With  the 
DG  guidance  law,  significant  improvements  to  the  missile  performance  in  the  rear 
quadrant,  tail  chase  scenario  were  achieved. 

For  the  non-maneuvering  target  scenarios,  it  was  shown  that  while  all  three 
guidance  laws  performed  similarly  under  no  noise  conditions,  when  noise  was  added  to 
the  system,  both  APN  and  DG  suffered  greater  perfonnance  degradation  as  compared  to 
PN.  This  is  intuitive  as  the  more  complex  guidance  laws  are  expected  to  be  affected  by 
noise  to  a  larger  extent. 

On  the  contrary,  when  noise  was  injected  for  the  maneuvering  scenarios,  the  DG 
law  was  shown  to  be  the  most  robust  performer  instead  of  PN,  while  APN  suffered  the 
greatest  perfonnance  degradations,  especially  around  the  60°  heading  region.  The  DG 
guidance  maintained  its  perfonnance  advantage  even  with  the  addition  of  noise  and  is 
shown  to  be  an  overall  superior  perfonning  guidance  law  that  can  be  relied  upon  in  both 
maneuvering  and  non-maneuvering  scenarios. 
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B.  FUTURE  RESEARCH 

Encouraged  by  the  results  of  the  current  study,  future  research  can  address  the 
following: 

•  Increasing  model  fidelity  with  more  accurate  sensor  and  controller  models. 

•  Expanding  the  kinematic  boundary  concept  to  three  dimensional  scenarios 
and  considering  a  broader  range  of  engagement  scenarios. 

•  Comparing  other  guidance  laws  against  the  ones  considered  in  this  study. 

•  Having  a  more  thorough  study  in  noise  implementation  and  specific 
filtering  algorithms. 
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APPENDIX  A.  SIMULINK  MODELS 


The  SIMULINK1'  block  diagrams  used  in  the  simulation  are  presented  in  this 
appendix.  The  simulation  is  created  as  a  single  model  and  functions  that  are  not  required 
can  be  modified  or  switched  off  as  required. 
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Output:  [  nx; 

ny: 

nz]: 


SIMPLIFIED  6D0F  SIMULATION 

>MODEL 
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Flat  Earth  6D0F  Dynamics 
from  Stevens  &  Lewis 
"Aircraft  Control  &  Simulation" 

>MODEL/Missile_Dynamics 


45 


13 


> 


d>7 

J 

CEH 

F  B 


> 


Fen: 

flatearthdyn 


> 


MATLAB 

Function 


dynamics 


13 


d>i - ► 

T  B 


GZh - ► 

mass 


>MODEL/Missile_Dynamics/Flat 
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position 


Initial  Condition: 
MSL  INIT 


Integrator 


KX) 

velocities 


-KD 

angular_rates 


■»m 

quartemions 


x  dot 


Earth  6  DOF  Model 


Gain: 


PD  controller  required  to  model 
missile  weather  vaning 

>MODEL/Moment  Feedback 


KjD 

moments 
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[X.  y;  2] 


cx>- 


Missile  Posn/Vel 


[x_dot;  y_dot;  z_dot] 


A 


[POSN1 


Simulation  of  Missile  IMU  and  Air  Data  Computer 

>MODEL/IMU_&_Air_Data_Computer2 


<D 
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as  truth  for  target 
turn  model 


Simulator  Seeker  Head  All  Data  Available  to  Guidance  Law 

User  can  select  which  to  use 

>MODEL/Seeker 
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cry*- 

LOS  rate  /  LOS 


Vert  LOS  noise 


Fen: 

atan2(u(3)-u(8),u(1  )-u(7)) 


du/dt 


MATLAB 

Function 


Derivative  theta 


theta  Jos 
Fen: 

atan2((u(5)-u(9)),sqrt((u(1)-u(7))A2+(u(3)-u(8))A2)) 


Gain 
Gain: 

57.296 

To  convert  radians  to  degress 


Seeke 


Goto 


>MODEL/Seeker/Angles_&_Velocity 


Ini 


>  Seeker 
r_outputsJo_workspace 
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51 


Missile  Guidance  Simulator  Target  Model 
Six  Dimensional  State  Vector 

>MODEL/Target_Model 
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This  block  sets  the  target 
turn_rate  in  terms  of 
centripetal  acceleration 
in  m/sA2 

For  no  turn,  set  to  0 


>MODEL/Target_Model/Target_Turn 
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Q> 

x  mhat 


Fen: 

abgfilter 


MATLAB  Fen 


Input: 

FILTSAMP 


MATLAB 

JA- 

- 

Function 

Zero-Order 

Hold 


>MODEL/Tracking  Filter 
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(JD-+ 

Ini 


Vx 


Ax 


Vy 


Ay 


Vz 


Input: 
XI  noise 


- ►€> 


Band-Limited 
White  Noise 


Input: 
X2  noise 


Band-Limited 
White  Noise3 


Input: 
X3  noise 


Band-Limited 
White  Noised 


Input: 
Y1  noise 


Band-Limited 
White  Noise  1 


Input: 
Y2  noise 


Band-Limited 
White  Noise4 


Input: 
Y3  noise 


Band-Limited 
White  Noise7 


Input: 
Z1  noise 


Band-Limited 
White  Noise2 


Input: 
Z2  noise 


Band-Limited 
White  Noise5 


Input: 
Z3  noise 


Band-Limited 
White  Noise8 


it - ►€> 


it - ►€> 


it - ►€> 


- ►€> 


it - ►£> 


- ►€> 


it - ►€> 


- ►£> 


>MODEL/Noise  Block 


+GD 

Outl 
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Upper/Lower  limit  set  to 
+/-  (MASS*30*GRAVITY) 


Simplified  Aerodynamic  Force  Generator 
Converts  Commanded  Acceleration  to  Force 

>MODEL/Aerodynamic_Force_Generator 


56 


2 

DragForce  4 - 

To  Workspace 


Fen: 

draginduced 


MATLAB 

Function 

4 


lnduced_Drag_Model 


Fen: 
dragparasitic 


MATLAB 

Function 


P  a  ra  s  it  e_D  rag_M  od  e  I 


16 


Control  Forces 


4 


4 


HD 

State 


HD 

Thrust 


Drag  Force  Model  for 
Simplified  6DOF  Flight  Model 

>MODEL/Drag_Model 
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Fen: 

sqrt(u(1)A2+u(2)A2+u(3)A2)/9.8045 


CT> 

G-Force 


Total  acceleration  in  G's 


>MODEL/G's 
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> 


GForceOut 


To  Workspace 


Fen: 

sqrt(u(2)A2+u(4)A2+u(6)A2) 


To_Workspace 


>MODEL/Tgt/Msl  Velocities 


59 


>MODEL/Sim  Stopper 
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APPENDIX  B.  MATLAB  CODE 


This  appendix  contains  the  list  of  script  files  and  function  codes  that  were  used  in 
the  simulation.  Table5  summarizes  the  code  names  and  their  purpose. 


FILENAME 

PURPOSE 

A.  Simulation  run  script  files 

Kbouter3.m 

Runs  and  computes  the  kinematic  boundary  for  a  particular  guidance  law 

NoiseStudy.m 

Runs  the  noise  simulation  for  a  100  repetitions 

Single  run.m 

Runs  a  single  simulation  for  a  particular  scenario 

B.  Simulation  initialization  script  files 

Missile  data.m 

Loads  the  pre-defmed  missile  constants  for  the  AMRAAM  model 

Thesis  init.m 

Loads  the  simulation  initialization  constants  and  variables 

C.  Simulation  guidance  law  function  files 

chingfanlin.m 

Implements  APN  guidance  law 

diffgeo.m 

Implements  DG  guidance  law 

propnavpt.m 

Implements  PN  guidance  law 

D.  Simulation  function  files 

abgfilter.m 

Calls  the  alpha-beta-gamma  filter  model 

alphabeta.m 

Calculates  the  AOA 

cdvmach.m 

Polyfits  the  data  curve  for  Cdi 

draginduced.m 

Calculates  the  induced  drag  force 

dragparasitic.m 

Calculates  the  parasitic  drag  force 

dynamic3D.m 

Runs  3D  target  dynamics 

flatearthdyn.m 

Runs  the  flat  earth  6DOF  dynamics 

formdrag.m 

Calculates  the  form  drag  of  the  missile 

machvalt.m 

Calculates  the  speed  of  Mach  one  at  particular  altitude 

q2euler.m 

Calculates  the  euler  angles  from  quarternions 

quarternion.m 

Calculates  the  quarternions  from  euler  angles 

quat2b.m 

Calculates  the  B  rotation  matrix  from  quarternions 

rhovalt.m 

Calculates  the  atmospheric  density  at  particular  altitude 

switchlimit.m 

Triggers  target  turn  rate  and  prevents  switch  back  to  zero 

tgo.m 

Calculates  the  time  to  go  for  the  missile  to  intercept  target 

thebigstop.m 

Stops  the  simulation  run  when  conditions  are  met 

Table  5.  Summary  of  MATLAB®  files  and  functions. 
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A. 


SIMULATION  RUN  SCRIPT  FILES 


o, 

o 

o. 


o 

o, 

o 

File : 

Kbouter3 .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Automatically  computes  a  kinematic  boundary  using  6 

o, 

o 

DOF  simulator  with  tracking  filter. 

o, 

o 

-  Streamlined  search  loops 

o, 

o 

-  Status  indicator 

o, 

o 

-  Saves  most  recent  data  to  disk 

o, 

o 

Inputs : 

none 

o, 

o 

Outputs : 

one  figure  of  kinematic  boundary 

o, 

o 

Process : 

streamlined  brute  force  search  algorithm 

o, 

o 

Assumptions : 

0, 

o 

CL 

Comments : 

o 

o, _ 

o 

clear 

clc 


% -  define  globals  - 

global  SWITCHFLAG  TURNFLAG 

% -  define  constants  - 

Thesis_init; 

STOPFLAG  =  0;  %  (1)  enable  display  of  simulation  stop  conditions 

%  (0)  disable  display  of  sim  stop  conditions 


% -  define  input  vector  - 

%  initialize  noise  variables 

%  Set  factor  to  0  if  noiseless  simulation  is  required,  else  set  to  any 


%  positive  integer  to  specify  noise 
%  mulitiplied  directly  into  the  powe 
%  white  noise  block 
factor  =  1; 

r  no ise  =  l*f actor;  % 

rdot_noise  =  le-2*factor;  % 

theta_noise  =  le-8*factor;  % 

thetadot_noise  =  le-8*factor;  % 

phi_noise  =  0*factor;  % 

phidot_noise  =  0*factor;  % 


.evel  desired.  Factor  will  be 
■  spectral  density  values  of  the 

Range  Noise 
Range  Rate  Noise 
Horizontal  LOS  angle 
Horizontal  LOS  angle  rate 
Vertical  LOS  angle 
Vertical  LOS  angle  rate 


%  set  noise  for  the  abg  filter 
abgfactor  =  factor; 

Xl_noise  =  l*abgfactor; 
X2_noise  =  le-2*abgfactor; 
X3_noise  =  2e-2*abgfactor; 
Yl_noise  =  l*abgfactor; 

Y2  noise  =  le-2*abgfactor; 
Y3_noise  =  2e-2*abgfactor; 


%  Target  X  posn  noise 
%  Target  X  velocity  noise 
%  Target  X  accel  noise 
%  Target  Y  posn  noise 
%  Target  Y  velocity  noise 
%  Target  Y  accel  noise 
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Zl_noise  =  O*abgfactor;  %  Target  Z  posn  noise 

Z2_noise  =  O*abgfactor;  %  Target  Z  velocity  noise 

Z3_noise  =  O*abgfactor;  %  Target  Z  accel  noise 

% -  initialize  variables  - 

TARGET  TURN  =6  %  set  target  turn  rate,  in  g's, 

%  default  =  0,  max  allowable  =  9. 

MIN  RNG  =  10000;  %  set  min  engagement  range  (10000m  default) 

DEGSTEP  =5;  %  set  heading  increment 

START_TIME  =  tic; 

MaxHit  =  zeros ( 1 0 0 0 , 5 ) ;  %  Initializing  1000x5  matrix  to  hold  MaxHit  data 
load  CURRENT 


% -  functions  - 

%  Start  in  tail  chase  step  to  head  on  by  <DEGSTEP>  degree  increments 
for  HEADING  =  0 : DEGSTEP : 1 8 0 


tic 

plotcount  =  1; 
runplot  =  zeros (100, 5) 


Initialize  100x5  matrix  to  hold  runplot 
data 

Preload  MIN  RNG  column  with  any  value 
greater  than  5  so  as  to  allow  for  correct 
index  search  later 

Initialize  1x2  matrix  to  hold  rangemax 
%  data 

rangemax ( 1 , 1 )  =  MIN  RNG;%  rangemax  is  a  1X2  matrix  to  store  rangemax 

%  and  min  Range  to  go  info 


runplot (:,1)  =  10; 


rangemax  =  [0  0]; 


disp ( [ ' Heading  ' , num2str (HEADING)  ,  '  deg ' ] ) 


show  heading  counter 


%  compute  target  speed  components 
XSPD  =  TGT_SPD*cos ( HEADING*pi / 1 8 0 )  ; 

YSPD  =  TGT_SPD*sin (HEADING*pi/180)  ; 

%  first  range  loop  step  by  10  km 

for  TGT  RNG  =  rangemax ( 1 , 1 )  :  10000  :  150000 

disp (['***  1 , num2str (TGT  RNG) , ' ,  10km  step  size*** ' ] ) 

%  set  initial  target  state 

TGT^INIT  =  [TGT_RNG;XSPD; 0; YSPD; -ALT; 0] ; 

TURNFLAG  =  0; 

SWITCHFLAG  =  0; 

XLAST  =  [TGT  RNG;0;0;0;0;  %  store  data  for  abgfilter.m  use 

0 ; -ALT ; 0 ; 0 ] ; 


%  call  simulation 
sim ( 'MODEL' ) 


%  save  run  data 

disp ( [ ' >>>  ’ , num2str (min (RangeToGo) ) , ' ,  min  Range<<< ' ] ) 
runplot (plotcount, : ) = [min (RangeToGo) ,0,0,0, TGT_RNG] ; 

%  score  run 
if  (min (RangeToGo) >5) 
break 
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end 

plotcount  =  plotcount+1; 

end 

INDEX  =  find (runplot ( : , 1) <=5)  ; 
if (INDEX) 

rangemax ( 1 , 1 )  =  runplot (max ( INDEX)  ,  5 )  ;  %  Store  rangemax  data 
rangemax ( 1 , 2 )  =  runplot (max ( INDEX) , 1 ) ;  %  Store  min  range2go  data 

end 

runplot  =  zeros (100, 5) ; 
runplot (:,1)  =  10; 
plotcount  =  1; 

%  main  search  loop  1km  step  size 

for  TGT  RNG  =  rangemax ( 1 , 1 ) +1 0 0 0  :  1000  :  rangemax ( 1 , 1 ) +9000 

disp (['***  '  , num2str (TGT  RNG) , ' ,  1km  step  size*** ' ] ) 

%  set  initial  target  state 

TGT_INIT  =  [TGT_RNG;XSPD; 0; YSPD; -ALT; 0] ; 

TURNFLAG  =  0; 

SWITCHFLAG  =  0; 

XLAST  =  [TGT  RNG;0;0;0;0;  %  store  data  for  abgfilter.m  use 

0 ; -ALT ; 0 ; 0 ] ; 

%  call  simulation 
sim ( 'MODEL' ) 

%  save  run  data 

disp ( [ ' >>>  , num2str (min (RangeToGo) ) , ' ,  min  Range<<< ' ] ) 

runplot (plotcount, : ) = [min (RangeToGo) ,0,0,0, TGT_RNG] ; 

%  score  run 
if  (min (RangeToGo) >5) 
break 

end 

plotcount  =  plotcount+1; 

end 

INDEX  =  find (runplot (:, 1) <=5) ; 
if (INDEX) 

rangemax ( 1 , 1 )  =  runplot (max ( INDEX) , 5 ) ;  %  Store  rangemax  data 
rangemax ( 1 , 2 )  =  runplot (max ( INDEX) , 1 ) ;  %  Store  min  range2go  data 

end 

runplot  =  zeros(100,5); 
runplot (:,1)  =  10; 
plotcount  =  1; 

%  main  search  loop  100m 

for  TGT  RNG  =  rangemax ( 1 , 1 ) +1 0 0  :  100  :  rangemax ( 1 , 1 ) +900 
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*  *  *  I 


disp (['***  ' , num2str (TGT  RNG) , ' ,  100m  step  size 
%set  initial  target  state 

TGT^INIT  =  [TGT_RNG;XSPD; 0; YSPD; -ALT; 0] ; 

TURNFLAG  =  0; 

SWITCHFLAG  =  0; 

XLAST  =  [TGT  RNG;0;0;0;0;  %  store  data  for  abgfilter.m  use 

0 ; -ALT ; 0 ; 0 ] ; 

%  call  simulation 
sim ( 'MODEL' ) 

%  save  run  data 

disp ( [ ' >>>  , num2str (min (RangeToGo) )  ,  '  ,  min  Range<<< ' ] ) 

runplot (plotcount, : ) = [min (RangeToGo) ,0,0,0, TGT_RNG] ; 

%  score  run 
if  (min (RangeToGo) >5) 
break 

end 

plotcount  =  plotcount+1; 

end 

INDEX  =  find (runplot (:, 1) <=5) ; 
if (INDEX) 

rangemax ( 1 , 1 )  =  runplot (max ( INDEX) , 5 ) ;  %  Store  rangemax  data 
rangemax ( 1 , 2 )  =  runplot (max ( INDEX) , 1 ) ;  %  Store  min  range2go  data 

end 

runplot  =  zeros (100, 5) ; 
runplot (:,1)  =  10; 
plotcount  =  1; 

%  main  search  loop  10m.  Note,  now  it  is  computing  the  full  output 
%  vector  for  each  run.  Starts  calculation  at  rangemax  so  as  to 
%  determine  full  output  vector  for  previous  rangemax. 
for  TGT  RNG  =  rangemax ( 1 , 1 )  :  10  : rangemax ( 1 , 1 ) +90 

disp(['***  ' , num2str (TGT  RNG),',  10m  step  size***']) 

%set  initial  target  state 

TGT^INIT  =  [TGT_RNG;XSPD; 0; YSPD; -ALT; 0] ; 

TURNFLAG  =  0; 

SWITCHFLAG  =  0; 

XLAST  =  [TGT  RNG;0;0;0;0;  %  store  data  for  abgfilter.m  use 

0 ; -ALT ; 0 ; 0 ] ; 

%  call  simulation 
sim ( 'MODEL' ) 

%  analyze  data  from  current  run 
TOUT  =  MissileOut  ( : , 14) ; 

INDX  =  find (RangeToGo==min (RangeToGo) ) ; 

IP  =  TOUT (INDX (1) ) ; 
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%  compute  cost  function  J=20*e (tf ) A2+integ (uA2 ) /200  and 
%  missile  divert 

u2  =  (OmegaOut ( : , 1 )  . A2+0mega0ut ( : , 2 )  . A2 )  ; 
integral  =  0; 
for  ii=2 : INDX 

integral  =  integral+ (TOUT (ii) -TOUT (ii-1) ) *u2 (ii-1) ; 

end 

J  =  20*min (RangeToGo) A2+integral/1000; 

%  save  run  data  [miss  dist,  cost,  divert,  time,  max  range] 
disp ( [ ' >>>  , num2str (min (RangeToGo) )  ,  ' ,  min  Range<<< ' ] ) 

runplot (plotcount,  : )  =  [min (RangeToGo) , J, integral,  IP,  TGT  RNG]; 

if  (min (RangeToGo) >5) 
break 

end 

plotcount  =  plotcount+1; 

end 

INDEX  =  find (runplot (:, 1) <=5) ; 
if (INDEX) 

rangemax ( 1 , 1 )  =  runplot (max ( INDEX) , 5 ) ;  %  Store  rangemax  data 
rangemax ( 1 , 2 )  =  runplot (max ( INDEX) , 1 ) ;  %  Store  min  range2go  data 

end 

if  ( isempty ( INDEX) ) 

MaxHit (HEADING+1 , : )  =  [rangemax (1, 2) , 0, 0, 0, rangemax (1, 1) ] ; 

else 

MaxHit (HEADING+1, : )  =  runplot (max ( INDEX) ,:) ; 

end 

%  save  data  to  disk 
save  CURRENT  MaxHit 
toe 

%  note  for  some  guidance  laws,  the  down  step  here 

%  must  be  2  or  more - I 

MIN  RNG  =  10000* (floor (rangemax (1, 1) /10000) -1) ; 
if  (MIN_RNG  <=  10000) 

MIN_RNG  =  10000; 

end 

%  MIN_RNG  =  10000; 

end 

toe (START_TIME) 

END  TIME  =  toe (START  TIME); 


O,  O, 
o  o 

%  plot  the  graph 

%  0  deg  represents  tail  chase  scenario 

%  180  deg  represents  head  on  scenario 

rhol  =  MaxHit ( 1 : DEGSTEP  :  18 1 , 5 )  ; 

rhol  =  [rhol ; flipud (rhol )] ; 

theta  =  0 :DEGSTEP: 180; 

theta  =  pi/180*theta; 
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theta  =  [theta, -l*fliplr (theta) ] ' 
figure ( 5 ) 

polar  (theta, rhol) 


o 

o, 

o 

File : 

Noise  Study. m 

o, 

o 

Name : 

CPT  Daniel  Perh 

g, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Automatically  runs  100  simulation  runs  at  a 

particular 

o, 

o 

point  with  random  noise  input  to  determine  noise 

o, 

o 

effects  on  guidance  laws. 

o, 

o 

Inputs : 

none 

g, 

o 

Outputs : 

plot  of  abg  filter  data  on  TGT  posn  for  the  last  run 

o, 

o 

Process : 

0, 

o 

Assumptions : 

o, 

o 

Comments : 

NOTE:  Must  switch  simulation  solver  to  a  fixed  step 

o, 

o 

o, 

o 

o. 

solver  before  running  this  script  file. 

clear 

clc 


% -  define  globals  - 

global  SWITCHFLAG  TURNFLAG 

% -  define  constants  - 

%  initialize  simulation 
Thesis_init; 

STOPFLAG  =  0;  %  (1)  enable  display  of  simulation  stop  conditions 

%  (0)  disable  display  of  sim  stop  conditions 

%  initialize  variables  to  hold  data 
holdrange=zeros (1000, 1)  ; 
holdpos=zeros (1000, 3)  ; 

% -  define  input  vector  - 

%  initialize  noise  variables 

%  Set  factor  to  any  positive  integer  to  specify  noise  level  desired. 

%  Factor  will  be  mulitiplied  directly  into  the  power  spectral  density 
%  values  of  the  white  noise  block 


factor  =  4.7 
r  noise  =  l*f actor ; 

O, 

o 

Range  Noise 

rdot  noise  =  le-2*factor; 

o, 

o 

Range  Rate  Noise 

theta  noise  =  le- 8* factor; 

Q, 

O 

Horizontal  LOS  angle 

thetadot  noise  =  le-8*factor; 

q, 

o 

Horizontal  LOS  angle  rate 

phi  noise  =  0*factor; 

o, 

o 

Vertical  LOS  angle 

phidot  noise  =  0*factor; 

o, 

o 

Vertical  LOS  angle  rate 

%  set  noise  for  the  abg  filter 
abgfactor  =  factor; 

XI  noise  =  l*abgfactor; 

g, 

o 

Target  X  posn  noise 

X2  noise  =  le-2*abgfactor; 

g, 

o 

Target  X  velocity  noise 

X3  noise  =  2e-2*abgfactor; 

g, 

o 

Target  X  accel  noise 

Y1  noise  =  l*abgfactor; 

g, 

o 

Target  Y  posn  noise 
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Y2_noise  =  le-2*abgfactor; 

Y3_noise  =  2e-2*abgfactor; 

Zl_noise  =  O*abgfactor; 

Z2_noise  =  O*abgfactor; 

Z3_noise  =  O*abgfactor; 

% -  initialize  variables 

%  initialize  target 
TGT_RNG  =  79900; 

HEADING  =  135; 

TARGE T_TURN  =  6; 

XSPD  =  TGT_SPD*cos (HEADING*pi/180)  ; 

YSPD  =  TGT_SPD*sin (HEADING*pi/180)  ; 

tic 

% -  functions  - 

%  100  realizations 
for  numloops=l : 100 
disp (numloops) 

TGT^INIT  =  [TGT_RNG; XSPD; 0;YSPD; -ALT; 0] ; 

TURNFLAG  =  0; 

SWITCHFLAG  =  0; 

XLAST  =  [TGT  RNG;0;0;0;0;  %  store  data  for  abgfilter.m  use.  Must 

0;-ALT;0;0];  %  be  reset  to  baseline  value  before 

%  each  simulation  call 

sim ( 'MODEL' ) 

%  analyze  data  from  current  run 
disp (min (RangeToGo) ) 

holdrange (numloops,  : ) =min (RangeToGo)  ; 
idx=find (RangeToGo==min (RangeToGo) )  ; 

holdpos (numloops , : ) =MissileOut ( idx, 1 : 3 ) -TgtOut ( idx, 1 : 2 : 5 ) ; 
save  NOISY  holdrange  holdpos 

end 

%  Calculate  and  display  mean  miss  distance  and  the  standard  deviation 
missdistance=mean (holdrange) 
sigmadistance=std (holdrange) 
toe 

%  Plot  target  position  based  on  abg  filter  results  for  the  last  run 
figure  (2 ) 

plot (XLAST_Data {:,!), XLAST_Data ( : , 4 ) ) 
title  (  ' Tgt  Posn  based  on  abg  filter') 
xlabel('X  Posn') 
ylabel('Y  Posn') 


%  Target  Y  velocity  noise 
%  Target  Y  accel  noise 
%  Target  Z  posn  noise 
%  Target  Z  velocity  noise 
%  Target  Z  accel  noise 
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o. 


o 

o, 

o 

File : 

Single  run.m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Runs  a  single  iteration  of  a  particular  engagement 

o, 

o 

scenario  with  user  specified  inputs 

o, 

o 

Inputs : 

Target  Parameters  (Alt,  Turn  rate,  Mach,  Hdg,  Rng) 

o, 

o 

Outputs : 

Min  Range  (min  miss  dist) 

o, 

o 

Graphs  of  (LOS  Range  to  Go  vs  Time) 

o, 

o 

(Engagement  Geometry) 

o, 

o 

(Msl  &  Tgt  Speed  vs  Time) 

0, 

o 

(Msl  &  Tgt  Accel  (g)  vs  Time) 

o, 

o 

(Guidance  command  output  vs  Time) 

o, 

o 

(Msl  Forces  vs  Time) 

o, 

o 

Process : 

o, 

o 

Assumptions : 

o, 

o 

CL 

Comments : 

o 

o,  _ 

o 

clear; 

clc; 


% -  define  globals  - 

global  STOPFLAG  SWITCHFLAG  TURNFLAG 

% -  define  constants  - 

Thesis_init; 

STOPFLAG  =  0;  %  (1)  enable  display  of  simulation  stop  conditions 

%  (0)  disable  display  of  sim  stop  conditions 

% -  define  input  vector  - 

%  initialize  noise  variables 

%  Set  factor  to  0  if  noiseless  simulation  is  required,  else  set  to  any 
%  positive  integer  to  specify  noise  level  desired.  Factor  will  be 
%  mulitiplied  directly  into  the  power  spectral  density  values  of  the 
%  white  noise  block 
factor  =  0 
r_noise  =  l*factor; 
rdot  noise  =  le-2*factor; 
theta_noise  =  le-8*factor; 
thetadot_noise  =  le-8*factor; 
phi  noise  =  0*factor; 
phidot_noise  =  0*factor; 

%  set  noise  for  the  abg  filter 
abgfactor  =  factor; 

Xl_noise  =  l*abgfactor; 

X2_noise  =  le-2*abgfactor; 

X3_noise  =  2e-2*abgfactor; 

Yl_noise  =  l*abgfactor; 

Y2_noise  =  le-2*abgfactor; 

Y3_noise  =  2e-2*abgfactor; 


Range  Noise 
Range  Rate  Noise 
Horizontal  LOS  angle 
Horizontal  LOS  angle  rate 
Vertical  LOS  angle 
Vertical  LOS  angle  rate 


Target  X  posn  noise 
Target  X  velocity  noise 
Target  X  accel  noise 
Target  Y  posn  noise 
Target  Y  velocity  noise 
Target  Y  accel  noise 
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Zl_noise  =  O*abgfactor;  %  Target  Z  posn  noise 

Z2_noise  =  O*abgfactor;  %  Target  Z  velocity  noise 

Z3_noise  =  O*abgfactor;  %  Target  Z  accel  noise 

% -  initialize  variables  - 

%  set  target  parameters 

ALT  =  6000;  %  default  co-altitude  in  metres 

TARGET  TURN  =6;  %  set  target  turn  rate,  in  g's, 

%  default  =  0,  max  allowable  =  9. 

%  set  target  speed 

TGT_MACH  =  0.83;  %  user  sets  Mach  #  for  target 

TGT_SPD  =  TGT  MACH*machvalt (ALT)  ;  %  machine  computes  speed 

%  define  scenario  variables 
TGT  HDG  =40;  %  heading  of 

%  heading  of 

TGT_RNG  =20000;  %  set  target 

% -  functions  - 

tic 

disp ([' Heading  =  ' , num2str (TGT  HDG),'  degrees']) 
disp ([' Target  Range  =  ' , num2str (TGT  RNG/1000),'  km']) 

%  compute  target  speed  components 
XSPD  =  TGT_SPD*cos (TGT_HDG*pi/180)  ; 

YSPD  =  TGT_SPD*sin (TGT_HDG*pi/180)  ; 

TGT^INIT  =  [TGT_RNG;XSPD; 0; YSPD; -ALT; 0] ; 

TURNFLAG  =  0; 

SWITCHFLAG  =  0; 

XLAST  =  [TGT  RNG;0;0;0;0;  %  store  data  for  abgfilter.m  use 

0 ; -ALT ; 0 ; 0 ] ; 

sim ( ' MODEL2 ' ) 
toe 

disp (['Min  Range  ' , num2str (min (RangeToGo) )  ,  '  m']) 


0  represents  tail  chase, 

180  represents  head  on  geometry 
range 


o,  o, 
o  o 

% - plot  graphs - 

time  =  rem(now,l); 

hr  =  f loor (time*24 ) ; 

mins  =  f loor (rem ( time*24 , 1 ) * 60 )  ; 

timestr  =  [ '  ' , num2str (hr) , ' : ' , num2str (mins) ] ; 


%  Missile  to  target  distance  (Range-to-go  vs  Time) 
range  =  RangeToGo; 
t  =  MissileOut ( : , 14 ) ; 
t  disc  =  0 : FILTSAMP :max (t) ; 


71 


index  =  find (range==min (range) ) ; 
ip  =  t ( index ( 1 ) ) ; 

figure  ( 1 ) 
subplot (2 , 1,2) 
plot (t, range) 

title ('LOS  Range  to  Go  vs  Time') 

xlabel ( [ ' time  (seconds)  , date, timestr] ) 

ylabel('LOS  Range  (meters) ') 

%  engagement  geometry 
subplot  (2 , 1,1) 

plot (TgtOut ( : , 1 ) , TgtOut ( : , 3 ) ,  ' :  '  ,  MissileOut ( : , 1 ) , MissileOut ( : , 2 ) ) 
axis  equal 

outtextl  =  ['time:  , num2str (ip) , '  seconds']; 
outtext2  =  ['range:  ', num2str (min (range) ), '  meters' ] ; 
text (100, 10000, ' Intercept  at:') 
text (100, 7000, outtextl) 
text (100, 4000, out text 2) 

title (' Engagement  Geometry') 
xlabel ('x  (meters) ') 
ylabel ( ' y  (meters) ' ) 

legend ( ' Target ' , ' Missile ' , 'Location ' , ' Best ' ) 


o,  o, 
o  o 

%  Missile  and  Target  Velocities 

Target_Spd  =  sqrt (TgtOut ( :  ,  2)  . A2+TgtOut ( :  ,  4)  . A2+TgtOut  (  : , 6)  .  A2)  ; 
figure (2 ) 

plot (t, Missile Speed, ' b- ' , t, Targe t_Spd, ' k : ' ) 
title (' Missile/Target  Speed  vs  Time') 
ylabel ( ' Speed  (m/s) ' ) 
xlabel ('Time  (s) ') 

legend ( ' Missile  1 ,  ' Target ' ,  'Location ' ,  ' Best ' ) 


o,  o. 
o  o 

%  Missile  &  Target  Accelerations  in  g's 
gforce  =  sqrt (AccelOut ( : , 1 ) . A2+AccelOut ( : , 2 ) . A2  ... 
+AccelOut (:,3) .A2) . / 9 . 8045 ; 

figure  ( 3 ) 
subplot  (2,1,1) 
plot (t, gforce) 

title (' Missile  Accelerations  (g)  vs  Time') 

ylabel (' Missile  Acceleration  (g)  ') 
axis([0  round (max (t) )  0  50]) 

%  compute  cost  function  J=20*e (tf ) A2+integ (uA2 ) /200 
u2  =  (OmegaOut ( : , 1 ) . A2+OmegaOut ( : , 2 ) . A2 ) ; 
integral  =  0; 
for  ii  =  2 : index 

integral  =  integral! (t (ii) -t ( i i — 1 ) ) *u2 (ii-1) ; 

end 

J  =  20*min (range) A2+integral/1000; 


72 


outtxt  =  ['Time  (s)  /  Missile  divert:  , num2str (integral) ] ; 

xlabel (outtxt) 

Tgtgf orce  =  ( sqrt ( Targe tVelAccel ( : , 2 ) . A2+TargetVelAccel ( : , 4 ) .  A2t . 

TargetVelAccel ( : , 6)  . A2) ) /9 . 8045; 
subplot (2 , 1,2) 
plot (t, Tgtgforce) 

title (' Target  Accelerations  (g)  vs  Time') 
xlabel ('Time  (s) ') 

ylabel (' Target  Acceleration  (g) ' ) 


o,  o, 
o  o 

%  guidance  command 
figure  ( 4 ) 

plot (t,  OmegaOut ( : , 1 ) , t, OmegaOut ( : , 2 ) ,  '  :  ' ) 

title (' Guidance  law  command  output  vs  Time') 

outtxt  =  ['Cost  J:  ' , num2str ( J)  ,  '  time  (seconds)']; 

xlabel ( [outtxt,  '  ' , date, timestr] ) 

ylabel ( ' n_c  (m/secA2) ') 

axis([0  round (max (t) )  -200  50]) 

legend ( ' n_c  y ' , ' n_c  z ' , ' Location ' , ' Best ' ) 


o,  o, 
o  o 

%  Missile  Force  vs  Time 
figure  ( 5 ) 

plot (t, ForcesOut ( : , 1) ,  b- ’ , t, ForcesOut ( : , 2) , ' k: ' ) 
title (' Missile  Forces  Fx/Fy  vs  Time') 
ylabel ( ' Force ' ) 
xlabel ('Time  (s) ') 

legend ( ' Fx ' , ' Fy ' , ' Location ' , ' Best ' ) 
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B. 


SIMULATION  INITIALIZATION  SCRIPT  FILES 


%%  AMRAAM  %% 

O, _ 

O 

o, _ 


o 

o, 

o 

File : 

Missile  data.m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584 

(R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Missile  data  for  AMRAAM.  Establishes 

missile 

o, 

o 

dimensions  for  use 

in  computing  aerodynamic  forces 

o, 

o 

and  moments .  Except 

where  noted,  all 

dimensions  in 

o, 

o 

MRS  system. 

0, 

o 

Inputs : 

o, 

o 

Outputs : 

various 

o, 

o 

Process : 

o, 

o 

Assumptions : 

o, 

o 

o. 

Comments : 

o 

o, 

o 


%  Missile  Name:  PSEUDO  AMRAAM 

% -  define  globals  - 

global  MASS  DIAM  LENGTH  XCG  XCPN  XCPB  XHL 

global  ST  SW  SPLAN  SREF  %  Areas 

global  JX  JY  JZ  %  Rotational  Inertia 


-  define  constants  - 

.ssile  body  dimensions 
156.8; 

0.1778; 

=  3.657; 

.8288; 

.  6769; 


% - mi 

MASS  = 

DIAM  = 

LENGTH 
XCG  =  1 
LN  =  0  . 

% -  missile  tailplane  dimensions 


mass,  may  be  time  varying 

diameter 

length 

initial  c.g.,  may  be  time  varying 
length  of  nose  cone 


XHL 

=  3.454; 

o, 

o 

hinge 

i  line  arm 

CRT 

=  0.4061; 

o, 

o 

tail 

root  chord 

CTT 

=  0.0676; 

o, 

o 

tail 

tip  chord 

TXT 

=  0.0676; 

o, 

o 

tail 

extension 

HT  = 

0.2286; 

o, 

o 

tail 

height 

O, 

O 

missile  wing 

dimensions  -- 

xw  = 

1.134; 

o, 

o 

wing 

to  radome  ■ 

CRW 

=  0.3554; 

o, 

o 

wing 

root  chord 

CTW 

=  0; 

o, 

o 

wing 

tip  chord 

WXT 

=  0; 

o, 

o 

wing 

extension 

HW  = 

0.1778; 

o, 

o 

wing 

height 

define  input  vector 
initialize  variables 
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% -  functions  - 

% -  Centers  of  pressure  - 

XCPN  =  0.67* LN ;  %  nose  CP 

XCPW  =  LN+XW+0 . 7*CRW-0 . 2*CTW;  %  wing  CP 

AN  =  0 . 67*LN*DIAM;  %  plan  area 

AB  =  (LENGTH-LN) *DIAM;  %  plan  area 

XCPB  =  (0 . 67*AN*LN+AB* (LN+0 . 5* (LENGTH-LN) ))/ (AN+AB) ;  %  body  CP 

% -  Area  computations  - 

SW  =  0 . 5*HW* (CTW+CRW) +CRW*WXT;  %  wing  area 

ST  =  0 . 5*HT* (CTT+CRT) +CRT*TXT;  %  tail  area 

SPLAN  =  (LENGTH-LN) *DIAM+0 . 67*LENGTH*DIAM;  %  body  and  nose  plan 

SREF  =  pi*DIAMA2/4;  %  missile  cross  sect 

area 

% -  Computing  the  inertial  matrix  - 

JX  =  MASS* ( (DIAM/2) A2) /2; 

JY  =  MASS* ( (LENGTHA2) /12+ ( (DIAM/2) A2) /4) +MASS* ( (LENGTH/2) -XCG) A2 
JZ  =  JY; 


of  nose 
of  body 


area 

onal 
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o. 


o 

q, 

o 

File : 

Thesis  init.m 

q, 

o 

Name : 

CPT  Daniel  Perh 

q, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

q, 

o 

32-bit  (win  32) 

q, 

o 

Date : 

08  July  2011 

q, 

o 

Description : 

This  script  file  initializes  thesis  work  missile 

q, 

o 

simulation 

q, 

o 

Inputs : 

q, 

o 

Outputs : 

q, 

o 

Process : 

q, 

o 

Assumptions : 

q, 

o 

0, 

Comments : 

o 

o, 

o 


% -  define  globals  - 

global  STOPFLAG  SATFLAG  SWITCHFLAG  TURNFLAG 


global  TARGET_TURN  XLAST  FILTSAMP 

% -  define  constants  - 

%  physical  constants 
OMEGA_X  =  7 . 292115e-5;  % 

GM_E  =  3.9860014el4;  % 

R_E  =  6.378164e6;  % 

F  =  1/298.257;  % 

OMEGA_E  =  [ OMEGA_X ; 0 ; 0 ] ;  % 

GRAVITY  =  9.8045;  % 


earth's  rotation  rate 

G*mass  of  earth 

radius  of  earth 

ellipsoidal  squash  factor 

earth's  rotational  velocity  vector 

gravitational  acceleration  constant 


% -  define  input  vector  - 

%  DEFAULT  INPUTS  REQUIRED  (Data  defined  here  can  be  superseded  by  the 
%  data  within  the  individual  simulation  script 

%  files) 


ALT  =  6000;  % 

TARGET_TURN  =0;  % 

O, 

o 

q, 

o 

q, 

o 

TGT_MACH  =0.83;  % 

TGT^SPD  =  TGT_MACH*machvalt (ALT)  ; 
MSL_MACH  =0.83;  % 

MSL_SPD  =  MSL_MACH*machvalt (ALT)  ; 
TGT_HDG  =70;  % 

q, 

o 

q, 

o 

TGT_RNG  =40000;  % 

q, 

o 


set  default  co-altitude  in  metres 
set  target  turn  rate,  in  g's, 
default  =  0,  max  allowable  =  9. 
Simulation  applies  tgt  turn  rate  8  s 
prior  to  impact 
user  sets  Mach  #  for  target 

user  sers  Mach  #  for  missile  at  launch 

target  heading  in  degrees.  0 
represents  tail  chase  and  180 
represents  head  on  geometry 
target  distance  from  initial  launch 
point 


%  REQUIRED  GLOBAL  VALUES 

STOPFLAG  =  0;  %  (1)  enable  display  of  simulation  stop 

%  conditions 

%  store  turn  rate  data,  always  set  to  (0) 

%  prevent  turn  rate  from  switching  back  to 
%  zero,  always  set  to  (0) 
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TURNFLAG  =  0; 
SWITCHFLAG  =  0; 


TMAX  =  200; 


set  simulation  max  run  time 


XLAST  =  [TGT_RNG; 0 ; 0 ; 0 ; 0 ;  %  store  data  for  abgfilter.m  use 

0 ; -ALT ; 0 ; 0 ] ; 

FILTSAMP  =0.01;  % 

o, 
o 

o, 
o 

o, 
o 

Q, 

O 

O, 

O 

O, 
o 

% -  initialize  variables  - 

%  missile  physical  parameters 
Missile_data; 

VB  =  [MSL  SPD;0;0];  %  Missile  initial  velocity  vector 

POSN  =  [0;0;-ALT];  %  initial  missile  position  vector,  note 

%  altitude  is  negative  in  NED  coord 

%  compute  Euler  angles  for  missile 

PSI  =  0*pi/180;  %  Varying  PSI  will  determine  the  angle  at 

%  which  the  missile  is  pointing  at  the 
%  target.  0  deg  will  be  pointing  straight 
%  at  the  target. 

THETA  =  0*pi/180; 

PHI  =  0*pi/180; 

% -  functions  - 

Q_0  =  quarternion (PHI, THETA, PSI) ; 

Q_0  =  Q_0/sqrt ( Q _ 0 (1) A2+Q_0 (2) A2+Q_0 (3) A2+Q_0 (4) A2) ; 

B  =  quat2b(Q_0); 

P  =  0*pi/180; 

Q  =  0*pi/180; 

R  =  0*pi/180; 

OMEGA_B  =  [ P; Q; R] ; 

%  missile  initial  state  vector 
MSL^INIT  =  [ POSN ; VB ; OMEGA_B ; Q_0 ] ; 

%  compute  target  speed  components 
XSPD  =  TGT_SPD*cos (TGT_HDG*pi/180)  ; 

YSPD  =  TGT_SPD*sin (TGT_HDG*pi/180)  ; 

%  target  initial  state  vector  [x; xdot ; y; y_dot ; z ; zdot ] 

TGT  INIT  =  [TGT_RNG;XSPD; 0; YSPD; -ALT; 0]; 


set  filter  sample  interval  for 
abgfliter.m  use.  Also  sets  the  Sampling 
time  for  the  white  noise  blocks  as  well 
as  the  step  size  interval  when  selecting 
fixed  step  size  simulation  solver.  Do  not 
set  larger  than  0.02  for  diffgeo  guidance 
law . 
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c. 


SIMULATION  GUIDANCE  LAW  FUNCTION  FILES 


function  [  y  ]  =  chingfanlin  (  u  ) 

%  CHINGFANLIN 

%  Computes  the  optimal  guidance  law  derived  by  Ching  Fan  Lin  pg  475  with 
%  drag  force  inputs  for  point  mass  simulation 


O, 

o 

o. 


o 

o, 

o 

File : 

chingfanlin .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

computes  APN  guidance 

law  from  Ching  Fan  Lin 

o, 

o 

Inputs  : 

seeker  output,  filter 

output,  accelerometer,  missile 

o, 

o 

timer 

o, 

o 

Outputs : 

command  accelerations. 

y  &  z  forces  for  drag 

o, 

o 

Process : 

o, 

o 

Assumptions : 

o, 

o 

O, 

Comments : 

o 

o, 

o 


% -  define  globals 

global  MASS  SATFLAG 


% -  define  constants 

Nprime  =  5; 

Nprimez  =  5; 


define  input  vector 


thetadot  = 

u  ( 1 )  ; 

phidot  = 

u  ( 2 )  ; 

los  = 

u  ( 3 )  ; 

philos  = 

u  ( 4  )  ; 

Vc  = 

-u  (5)  ; 

R  = 

u  ( 6 )  ; 

heading  = 

u  (7 )  ; 

Vm  = 

u  ( 8 )  ; 

Vmdot  = 

u  ( 9 )  ; 

phi  = 

u  ( 1 0 )  ; 

theta  = 

u(ll)  ; 

psi  = 

u  ( 1 2 )  ; 

tgt  state  = 

u  (13 : 21) 

time  = 

u  (22)  ; 

accel  in  = 

u  ( 2  3  :  2  5 )  ; 

initialize  variables 


functions 
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if  (Vc==0 ) 

tgo  =  le6; 

else 

tgo  =  R/Vc; 

end 

%  Compute  relative  state  estimate 
xhat  =  [R*cos(los); 

R*sin (los) ; 

R*sin (philos) ; 
tgt_state (2) -Vm*cos (psi) ; 
tgt_state (5) -Vm*sin (psi) ; 
tgt_state (8) -Vm*sin (theta)  ; 
tgt_state ( 3 ) ; 
tgt_state ( 6 ) ; 
tgt_state ( 9 ) ; 
accel_in ( 1 ) ; 
accel  in  (2 ) ; 
accel  in  ( 3 ) ]  ; 
if  time<2 . 0 

ny  =  Nprime*Vc* (thetadot) /cos (heading-los) ; 
nz  =  Nprimez*Vc* (phidot) -9 . 8045; 

else 

uc  =  (5/tgoA2) * [eye (3) , tgo*eye (3) ,  tgoA2/2*eye (3) ,  zeros (3) ] *xhat; 

ny  =  uc (2) ; 

nz  =  uc(3)  -  9.8045; 

end 

%  Define  acceleration  in  x-axis  as  0  as  it  will  be  determined  from  the 
%  thrust  and  drag  model  blocks 
nx  =  0  ; 

%  Saturation  of  forces  at  30  g's  will  be  done  in  the  Aerodynamic  Force 
%  Generator  Block 
y  =  [nx; ny; nz] ; 

end 
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function  [  y  ]  =  diffgeo  (  u  ) 

%  DIFFGEO 

%  Computes  the  differential  geometric  guidance  law  derived  by  Chaoyong 
Li 


O, 


o 

o, 

o 

o, 

o 

File : 

diffgeo .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

computes  Differential 

Geometry  guidance 

law  from 

o, 

o 

Chaoyong  Li 

o, 

o 

Inputs : 

seeker  output,  filter 

output,  accelerometer,  missile 

o, 

o 

timer 

o, 

o 

Outputs : 

command  accelerations. 

y  &  z  forces  for 

drag 

o, 

o 

Process : 

o, 

o 

Assumptions : 

o, 

o 

CL 

Comments : 

o 

O, 

% -  define  globals 

global  MASS  SATFLAG 


% -  define  constants  -- 

Nprime  =  5; 

Nprimez  =  5; 

% -  define  input  vector 


thetadot  = 

u  ( 1 )  ; 

phidot  = 

u  ( 2 )  ; 

los  = 

u  ( 3 )  ; 

philos  = 

u  ( 4  )  ; 

Vc  = 

-u  (5)  ; 

R  = 

u  ( 6 )  ; 

heading  = 

u  (7 )  ; 

Vm  = 

u  ( 8 )  ; 

Vmdot  = 

u  ( 9 )  ; 

phi  = 

u  ( 1 0 )  ; 

theta  = 

u(ll)  ; 

psi  = 

u ( 12 )  ; 

tgt  state  = 

u  ( 13  :  21 ) 

time  = 

u (22)  ; 

accel  in  = 

u  ( 2  3  :  2  5 )  ; 

% -  initialize  variables  - 

% -  functions  - 

tgt_accel  =  [u (15) ;u (18) ] ;  %  Target  acceleration  vector  in  X-Y 

%  plane 
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eta_t  =  atan2 (u (17) ,u (14) )  -  los;  %  Angle  between  velocity  vector  of 

%  target  and  the  los  vector  in 
%  horizontal  plane 

eta  tz  =  atan2 (u (2 1 ) , norm ( tgt  accel) ) ;  %  angle  eta  in  vertical  plane 
eta_m  =  psi  -  los;  %  Look  angle  of  the  missile 


%  Diff  Geometric  guidance  law 

ny  =  norm (tgt_accel) *cos (eta_t) /cos (eta_m)  +... 

Nprime*Vc* (thetadot) /cos (eta_m) ; 

%  vertical  acceleration  must  account  for  gravity 
nz  =  u (21) *cos (eta_tz) /cos (theta-philos)  +... 

Nprimez*Vc* (phidot) /cos (theta-philos)  -  9.8045; 

%  Define  acceleration  in  x-axis  as  0  as  it  will  be  determined  from  the 
%  thrust  and  drag  model  blocks 
nx  =  0  ; 

%  Saturation  of  forces  at  30  g's  will  be  done  in  the  Aerodynamic  Force 
%  Generator  Block 
y  =  [nx;  ny;  nz]  ; 

end 
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function  [  y  ]  =  propnavpt  (  u  ) 

%  PROPNAVPT 

%  Computes  the  exact  proportional  navigation  with  drag  force  inputs  for 
%  point  mass  simulation 


O, 

o 

o. 


"6 

%  File: 

propnavpt .m 

%  Name : 

CPT  Daniel  Perh 

%  Compiler: 

MatLab  v7. 11. 0.584  (R2010b) 

O, 

o 

32-bit  (win  32) 

%  Date : 

08  July  2011 

%  Description: 

Prop  nav  guidance  law  for  6  DOF 

flight  model. 

Computes 

o, 

o 

the  applied  forces  for  use  by  induced  drag  model. 

o, 

o 

Required  to  eliminate  algebraic 

loops  in  the 

simulation 

%  Inputs: 

[seeker  data,  IMU  data,  timer] 

%  Outputs: 

[command  accelerations,  applied 

forces ] 

%  Process: 

%  Assumptions: 

%  Comments: 

O, 

o 

o, 

o 


% -  define  globals 

global  MASS  SATFLAG 


define  constants 


Nprime  = 

5; 

Nprimez 

=  5; 

Q_ 

O 

define  input 

thetadot 

=  u  ( 1 )  ; 

phidot  = 

u  ( 2  )  ; 

los  = 

u  ( 3 )  ; 

philos  = 

u  ( 4  )  ; 

Vc  = 

-u (5)  ; 

theta  = 

u(ll)  ; 

psi  = 

u  (12)  ; 

% -  initialize  variables  - 

% -  functions  - 

%  classic  PN  guidance  law 

ny  =  Nprime*Vc* (thetadot) /cos (psi-los) ; 

%  vertical  acceleration  must  account  for  gravity 
nz  =  Nprimez*Vc* (phidot) /cos (theta-philos) -9 . 8045; 

%  Define  acceleration  in  x-axis  as  0  as  it  will  be  determined  from  the 
%  thrust  and  drag  model  blocks 
nx  =  0 ; 
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o\°  o\° 


Saturation  of  forces  at  30  g's  will  be  done  in  the  Aerodynamic  Force 
Generator  Block 
y  =  [nx; ny;  nz]  ; 

end 
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D. 


SIMULATION  FUNCTION  FILES 


function  [  y  ]  =  abgfilter  (  u  ) 

%  ABGFILTER 

%  Implements  an  alpha-beta-gamma  filter  as  outlined  in  Bar-Shalom  &  Li 
%  "Estimation  and  Tracking" 


O, 


o 

0, 

o 

o, 

o 

File : 

abgfilter .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Implements  a  9-dimensional  state  vector 

o, 

o 

alpha-beta-gamma  tracking  filter  for  use 

with 

o, 

o 

missile  guidance  laws  requiring  tracking 

filters 

o, 

o 

(Note:  Uses  global  XLAST  to  preserve  state  between 

o, 

o 

iterations ) 

o, 

o 

Inputs : 

measurements  (los, los  dot,R,R  dot). 

o, 

o 

missile  posn  (x,y,z) 

o, 

o 

Outputs : 

9-dimensional  estimate  of  target  state 

o, 

o 

(x, vx, ax, y, vy, ay, z , vz , az ) 

o, 

o 

Assumptions : 

o, 

o 

Comments : 

May  require  up  to  20  samples  to  stabilize 

from 

o, 

o 

0, 

o 

O, 

initialization 

% -  define  globals 

global  FILTSAMP  XLAST 


% -  define  constants  - 

% -  define  input  vector  - 

losdot  =  u ( 1 )  ; 
phidot  =  u (2 ) ; 

1  o  s  =  u  ( 3 )  ; 
phi  =  u (4) ; 
rdot  =  u  ( 5 ) ; 

R  =  u  ( 6 )  ; 
xm  =  u ( 7 ) ; 
ym  =  u  ( 8 )  ; 
zm  =  u  (  9 ) ; 

% -  initialize  variables  - 

%  compute  target  cartesian  coordinates 
xt  =  R*cos ( los ) +xm; 
yt  =  R*sin (los) +ym; 
zt  =  R*sin (phi) +zm; 

z  =  [xt; yt; zt] ; 
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%  set  noise  parameters 
sigmav  =  1; 
sigmaw  =  1; 

lamda  =  sigmav* (FILTSAMPA2 ) /sigmaw; 

%  set  filter  parameters  from  Bar-Shalom  &  Li  (Assumed  Numbers) 
f alpha  =  .9; 
fbeta  =  .9; 
f gamma  =  . 9 ; 

%  filter  matrices 

F  =  [1  FILTSAMP  FILTSAMPA2/2  zeros (1,6); 

0  1  FILTSAMP  zeros (1,6); 

001  zeros (1,6) ; 

0001  FILTSAMP  FILTSAMPA2/2  zeros (1,3); 

00001  FILTSAMP  zeros (1,3); 

000001  zeros (1,3); 

0000001  FILTSAMP  FILTSAMPA2 /2 ; 

00000001  FILTSAMP; 

000000001]; 

H  =  [100000000; 

000100000; 

000000100]; 

%  compute  steady  state  gains 

W  =  [falpha; fbeta/FILTSAMP; fgamma/ (2*FILTSAMPA2) ] ; 

%  build  gain  matrix 
P  =  [W  zeros (3,2) ; 

zeros (3,1)  W  zeros (3,1); 
zeros (3,2)  W] ; 

% -  functions  - 

%  run  filter 

xhat  =  F*XLAST; 

xhatl  =  xhat  +  P* ( z-H*xhat) ; 

XLAST  =  xhatl; 

y  =  xhatl; 

end 
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o\°  o\° 


function  [  y  ]  =  alphabeta (  u  ) 

ALPHABETA 

Computes  angles  of  attack  in  both  vertical  and  horizontal  planes 


O, 

o 

o. 


o 

o, 

o 

File : 

alphabeta .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

o, 

o 

Description : 

Computes  angles  of  attack  using  ATAN  formulation  in 
Bryson  "Control  of  Spacecraft  and  Aircraft" 

o, 

o 

Inputs : 

Missile  state 

o, 

o 

Outputs : 

Angles  of  attack  (alpha,  beta) 

o, 

o 

o, 

o 

Process : 
Assumptions : 

ATAN  formulation  of  Bryson 

o, 

o 

CL 

Comments : 

o 

o, 

o 


% -  define  globals  - 

% -  define  constants  - 

% -  define  input  vector  - 

v  =  [  u  ( 4 )  ;  u  ( 5 )  ;  u  ( 6 )  ]  ; 

% -  initialize  variables  - 

% -  functions  - 

%  (Equations  developed  in  Bryson) 

%  using  betal  for  sideslip  angle  to  avoid  problems  with  built  in  matlab 
%  function  'beta' 

alpha  =  atan2 (v (3) , sqrt (v (1) A2+v  (2) A2) )  ; 
betal  =  atan2 (v (2 ) , v (1 ) ) ; 

y  =  [alpha;betal] ; 

end 
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function  [  y  ]  =  cdvmach  (mach,  boost) 

%  CDVMACH 

%  Computes  approximation  of  zero  lift  drag  coefficient  vs.  mach  number 


File : 

cdvmach .m 

Name : 

CPT  Daniel  Perh 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

32-bit  (win  32) 

Date : 

08  July  2011 

Description : 

Computes  polynomial  fit  for  cdO  vs  Mach  number 

Inputs : 

mach  #  and  boost  status 

Outputs : 

cdO 

Process : 

Fit  on  data  from  Hutchins  EC4330  notes 

Assumptions : 

Comments : 

define  globals  - 
define  constants 


NoBoost  =  [-0.0014  0.0299  -0.2110  0.6256]; 
Boost  =  [-0.0012  0.0243  -0.1521  0.4044]; 

% -  define  input  vector  - 

% -  initialize  variables  - 

% -  functions  - 

if  (boost  &  (mach<0.7)) 
y=0 .15; 

end 

if  (-boost  &  (mach<0.7)) 
y=0 .25; 

end 

if  (boost  &  (mach>=0.7)  &  (mach<1.2)) 

y= (mach-0 . 7 ) *0 . 2  +  0.15; 

end 

if  (-boost  &  (mach>=0.7)  &  (mach<1.2)) 

y= (mach-0 . 7 ) *0 . 3  +  0.25; 

end 

if  ((mach>=1.2)  &  (boost~=0)) 

y=polyval (Boost,  mach) ; 

end 

if  ((mach>=1.2)  &  (boost==0)) 

y=polyval (NoBoost,  mach) ; 

end 

if  ( (mach>5  &  boost) ) 
y=0 . 10 ; 

end 

if  ((mach>6.4)  &  -boost) 

y=0 .132; 

end 

end 
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draginduced (  u  ) 


function  [  y  ]  = 

%  DRAGINDUCED 
%  Computes  the  induced  aerodynamic  drag  force 


o, 

o 

o. 


o 

o, 

o 

File : 

draginduced . m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

computes  induced  drag  for 

simplified  6DOF 

o, 

o 

Inputs : 

force  output  of  guidance 

law,  state 

o, 

o 

Outputs : 

drag  force 

o, 

o 

Process : 

work  backwards  to  CN  from 

forces 

o, 

o 

Assumptions : 

o, 

o 

O, 

Comments : 

o 

o, 

o 


% -  define  globals  - 

global  MASS  SREF 

% -  define  constants  - 

eAR  =1.5;  %  elliptical  eff  &  AR 


% -  define  input  vector  - 

Fy  =  u  (2) ; 

Fz  =  u (3) ; 

v2  =  u (7) A2+u (8) A2+u (9) A2;  % 

a  1 1  =  u  ( 6 )  ;  % 

% -  initialize  variables  - 

rho  =  rhovalt (abs (alt) ) ;  % 

mach  =  sqrt (v2 ) /machvalt (alt) ; 

Q  =  rho*v2/2;  % 

% -  functions  - 

if  (Q==0 ) 

Cny  =  0 ; 

Cnz  =  0; 

else 

Cny  =  Fy/ (Q*SREF) ;  % 

Cnz  =  Fz/ (Q*SREF) ;  % 

end 

Cdi  =  (CnyA2+Cnz A2 ) / (pi*eAR) ;  % 

if  (mach<l) 

Cdi  =  0 .25*sqrt (FyA2+FzA2) / (MASS* 

end 

y  =  Cdi*Q*SREF;  % 

end 


missile  velocity 
missile  alt 

atmospheric  density 
dynamic  pressure 


y  normal  coefficient 
z  normal  coefficient 

induced  drag  coefficient 

.8045);  %  subsonic  drag  equal  to 

%  max  CdO*applied  G  force 

drag  force 
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function  [  y  ]  =  dragparasitic  (  u  ) 

%  DRAGPARASITIC 

%  Computes  parasitic  aerodynamic  drag  force 


o, 

o 

o. 


o 

o, 

o 

File : 

dragparasitic .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Computes  parasitic  drag  after  breaking  apart  state 

o, 

o 

vector 

o, 

o 

Inputs : 

state  vector,  boost  status 

o, 

o 

Outputs : 

parasitic  drag  force 

o, 

o 

Process : 

o, 

o 

Assumptions : 

o, 

o 

CL 

Comments : 

o 

o, 

o 


% -  define  globals  - 

global  SREF 

% -  define  constants  - 

% - define  input  vector - 

vel2  =  u  (4) A2+u (5) A2+u (6) A2; 
alt  =  u ( 3 ) ; 

boost  =  u ( 14 ) ; 

% -  initialize  variables  - 

% -  functions  - 

y  =  formdrag (SREF,  alt,  vel2,  boost); 

end 
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function  [  y  ]  =  dynamic3D  (  u  ) 

%  DYNAMIC3D 

%  Computes  the  motion  dynamics  for  a  body  (target)  in  three  dimensions 


O, 

o 

o. 


o 

o, 

o 

File : 

dynamic3d .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7 . 11 . 0 

.584  (R2010b) 

o, 

o 

32-bit  (win  32 

) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

target  motion 

dynamics 

o, 

o 

Inputs : 

target  state. 

turn  rate  input 

o, 

o 

Outputs : 

derivative  of 

target  state 

o, 

o 

Process : 

o, 

o 

Assumptions : 

o, 

o 

o. 

Comments : 

o 

o, 

o 


o 

-  define 

globals  - 

o, 

o 

-  define 

constants  - 

0, 

o 

-  define 

input  vector  - 

Ac 

=  u  ( 1 )  ; 

%  Centripetal  Acceleration 

x  =  u  ( 2  )  ; 
xdot  =  u  ( 3 ) ; 
y  =  u  (4) ; 
ydot  =  u  ( 5 ) ; 
z  =  u  ( 6 )  ; 
zdot  =  u  ( 7 ) ; 

% -  initialize  variables  - 

% -  functions  - 

TgtSpd  =  sqrt (xdot A2+ydot A2+zdot A2 ) ; 

omega  =  Ac/TgtSpd;  %  where  Ac  =  TgtSpdA2  /  r 

%  &  omega  =  TgtSpd  /  r 

y  =  [xdot; 

-omega*ydot; 

ydot; 

omega*xdot; 

zdot; 

0]  ; 

end 
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function  [  y  ]  =  flatearthdyn  (  u  ) 

%  FLATEARTHDYN 

%  Computes  motion  dynamics  for  6  DOF  flat  earth  model 


O, 

o 

o. 


o 

o, 

o 

File : 

flatearthdyn . m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Computes  6  DOF  dynamics  for  flat  earth  using 

o, 

o 

quarternion  formulation 

o, 

o 

Inputs : 

o, 

o 

Outputs : 

derivative  of  state  vector 

o, 

o 

Process : 

Steven  &  Lewis 

o, 

o 

Assumptions : 

o, 

o 

O, 

Comments : 

o 

o, 

o 


define  globals  - 
define  constants 


% -  define  input  vector  - 

p  =  [  u  ( 1 )  ;  u  ( 2 )  ;  u  ( 3 )  ]  ; 
v_b  =  [  u  ( 4 )  ;  u  ( 5 )  ;  u  ( 6 )  ]  ; 
omega_b  =  [ u ( 7 ) ; u ( 8 ) ; u ( 9 ) ] ; 

P  =  u  ( 7  )  ;  Q  =  u  (  8  )  ;  R  =  u  (  9 )  ; 

q  =  [  u  ( 1 0  )  ;  u  ( 1 1 )  ;  u  ( 1 2  )  ;  u  ( 1 3 )  ]  ; 

magq  =  sqrt (q (1) A2+q (2) A2+q (3) A2+q (4) A2) ; 


q  =  q/magq; 

x  =  [p; v_b; omega_b; q] ; 

J  =  [u (14)  0  0; 

0  u  ( 1 5 )  0  ; 

0  0  u  ( 1 6 )  ]  ; 

F_B  =  [u  (17) ; u (18) ;u(19) ] ; 

T_B  =  [u (20 ) ; u ( 2 1 ) ; u ( 2  2 ) ] ; 

g  =  [0;  0;  9.8045] ; 

m  =  u  (23)  ; 

% -  initialize  variables 

%  Compute  rotation  matrices 
B  =  quat2b (q) ; 


%  inertial  matrix 

%  Forces 
%  Torques 

%  Not  using  external  gravity  model 
%  Mass 
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OMEGA  B 


OMEGA  B  = 

[0 

R 

o 

1 

Q; 

-P; 

-Q 

P 

0] 

r 

OMEGA  q  = 

[0 

P 

Q 

R; 

-P 

0 

-R 

Q; 

-Q 

R 

0 

-P; 

-R 

-Q 

P 

0] 

% -  functions 

y  =  [ zeros  ( 3 ) 
zeros  ( 3 ) 
zeros  ( 3 ) 
zeros (4,3) 


B' 

-OMEGA_B 
zeros ( 3 ) 
zeros (4,3) 


zeros ( 3 ) 
zeros ( 3 ) 

-l*inv ( J) *OMEGA_B* J 
zeros (4,3) 


zeros (3,4) ; 
zeros (3,4) ; 
zeros (3,4) ; 
(1/2) *OMEGA_q] 


y  =  y*x; 


y  =  y+ [ zeros  ( 3 , 1 )  ; 

B*g+ (1/m) *F_B; 
inv ( J) *T_B; 
zeros  (4,1)]; 


end 
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function  [  y  ]  =  formdrag  (A,  alt,  vel2,  boost  ) 

%  FORMDRAG 

%  Computes  form  drag  for  a  missile  with  frontal  area  A  in  a  standard 
%  atmosphere.  Uses  MACHVALT,  CDVMACH ,  RHOVALT 


O, 

o 

o. 


o 

o, 

o 

File : 

formdrag .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7 . 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Computs  the  form  drag  for  a  missile  with  frontal 

o, 

o 

area  A  in  a  standard  atmosphere 

o, 

o 

Inputs : 

area,  altitude.  VA2,  boost  on/off 

o, 

o 

Outputs : 

parasitic  drag  force 

o, 

o 

Process : 

o, 

o 

Assumptions : 

0, 

o 

o. 

Comments : 

o 

o, 

o 


% -  define  globals  - 

% -  define  constants  - 

% - define  input  vector - 

% -  initialize  variables  - 

rho  =  rhovalt (alt) ; 

mach  =  (vel2) A (1/2) /machvalt (alt) ; 

% -  functions  - 

if  (mach>100) 

mach  =  0.83; 

end 

Cd  =  cdvmach (mach, boost)  ; 
y  =  rho*vel2*Cd*A/2 ; 
end 


93 


function  [  y  ]  =  machvalt  (  alt  ) 

%  MACHVALT 

%  Computes  the  linear  approximation  for  a  given  altitude  in  meters/sec 
%  based  on  standard  ICAO  atmosphere 


o, 

o 

o. 


o 

o, 

o 

File : 

machvalt .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Computes  the  linear  approximation  to  Mach  1  for 

o, 

o 

standard  ICAO  atmosphere 

o, 

o 

Inputs : 

altitude 

o, 

o 

Outputs : 

Mach  1 

o, 

o 

Process : 

o, 

o 

Assumptions : 

0, 

o 

o. 

Comments : 

o 

o, 

o 


% -  define  globals  - 

% -  define  constants  - 

Machl  =  [-0.0041  340.3] ; 

Mach2  =  295.1; 

Mach3  =  [0.00067  281.7] ; 

% -  define  input  vector  - 

% -  initialize  variables  - 

alt  =  abs (alt) ;  %  account  for  NED  coords 

% -  functions  - 

if  (alt<11000 ) 

y  =  polyval (Machl , alt) ; 

else 

if  (alt>20000 ) 

y  =  polyval (Mach3, alt)  ; 

else 

y  =  Mach2; 

end 

end 

end 
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q2euler (  u  ) 


function  [  y  ]  = 

%  Q2EULER 

%  Computes  the  Euler  angles  from  quarternions 


O, 

o 

o. 


o 

o, 

o 

File : 

q2euler .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584 

(R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Computes  the  euler 

angles  from  the  quarternions 

o, 

o 

Inputs : 

Quarternions 

o, 

o 

Outputs : 

Euler  Angles 

o, 

o 

Process : 

Kuiper 

o, 

o 

Assumptions : 

o, 

o 

CL 

Comments : 

o 

o, 

o 


% -  define  globals  - 

% -  define  constants  - 

% - define  input  vector - 

qO  =  u  ( 1 )  ; 
ql  =  u  (2)  ; 
q2  =  u  ( 3 )  ; 
q3  =  u  ( 4 )  ; 

% -  initialize  variables  - 

% -  functions  - 

%  convert  quarternions  to  euler  angles 

mil  =  2*qOA2+2*qlA2-l; 

ml 2  =  2*ql*q2+2*q0*q3; 

ml 3  =  2*ql*q3-2*q0*q2; 

m23  =  2*q2*q3+2*q0*ql ; 

m3 3  =  2*qOA2+2*q3A2-l; 

psi  =  atan2 (ml2 , mil ) ; 
theta  =  asin(-ml3); 

%  correct  for  singularity  in  pitch 
if  (isreal (theta) ) 
theta  =  theta; 

else 

theta  =  sign (-ml3) *pi/2; 

end 

phi  =  atan2 (m23,m33) ; 
y  =  [phi,  theta,  psi] ; 

end 
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function  [  y  ]  =  quarternion  (  phi , theta, psi  ) 

%  QUARTERNION 

%  Computes  the  quarternions  from  Euler  angles  in  radians 


File : 

quarternion .m 

Name : 

CPT  Daniel  Perh 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

32-bit  (win  32) 

Date : 

08  July  2011 

Description : 

Computes  the  quarternions  from  euler  angles 

Inputs : 

Euler  angles  in  radians 

Outputs : 

Quarternions 

Process : 

Kuiper 

Assumptions : 

Comments : 

% -  define  globals  - 

% -  define  constants  - 

% -  define  input  vector  - 

% -  initialize  variables  - 

% -  functions  - 

%  Quarternion  equations 

qO  =  cos (phi/2 ) *cos (theta/2 ) *cos (psi/2 )..  . 

+  sin  (phi/2) *sin (theta/2) *sin (psi/2)  ; 
ql  =  sin (phi/2 ) *cos (theta/2 ) *cos (psi/2 ).. . 

-cos (phi/2) *sin (theta/2) *sin (psi/2)  ; 
q2  =  cos (phi/2 ) *sin (theta/2 ) *cos (psi/2  )..  . 

+  sin (phi/2) *cos (theta/2) *sin (psi/2)  ; 
q3  =  cos (phi/2 ) *cos (theta/2 ) *sin (psi/2 ).. . 

-sin (phi/2) *sin (theta/2) *cos (psi/2)  ; 

y  =  [qO;  ql;  q2 ;  q3]; 

end 
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function  [  y  ]  =  quat2b  (  u  ) 

%  QUAT2B 

%  Computes  rotation  matrix  from  quarternions 


File : 

quat2b .m 

Name : 

CPT  Daniel  Perh 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

32-bit  (win  32) 

Date : 

08  July  2011 

Description : 

Computes  rotation  matrix  from  quarternions 

Inputs : 

Quarternion 

Outputs : 

Rotation  Matrix  B 

Process : 

Kuiper 

Assumptions : 

Comments : 

% -  define  globals  - 

% -  define  constants  - 

% -  define  input  vector  - 

qO  =  u  ( 1 )  ; 
ql  =  u  (2)  ; 
q2  =  u  ( 3 )  ; 
q3  =  u  ( 4 )  ; 

% -  initialize  variables  - 

% -  functions  - 

2* (ql*q3-q0*q2 ) ; 

2* (q2*q3+q0*ql) ; 
q0A2-qlA2-q2A2+q3A2] ; 

end 


y  =  [q0A2+qlA2-q2A2-q3A2 
2* (ql*q2-q0*q3) 

2* (ql*q3+q0*q2 ) 


2* (ql*q2+q0*q3) 
qO  A2 -ql A2+q2  A2 -q3  A2 
2* (q2*q3-q0*ql ) 
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function  [  y  ]  =  rhovalt  (  alt  ) 

%  RHOVALT 

%  Computes  the  atmospheric  density  vs  altitude  for  ICAO  standard 
atmosphere 


O, 

o 

o. 


o 

o, 

o 

File : 

rhovalt .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Computes  the  atmospheric  density  from  ICAO  standard 

o, 

o 

atmosphere.  Exponential  model 

o, 

o 

Inputs : 

Altitude 

o, 

o 

Outputs : 

rho 

o, 

o 

Process : 

o, 

o 

Assumptions : 

0, 

o 

o. 

Comments : 

o 

o, 

o 


% -  define  globals  - 

% -  define  constants  - 

% -  define  input  vector  - 

% -  initialize  variables  - 

alt  =  abs (alt) ;  %  account  for  NED  coord 

% -  functions  - 

if  alt>9144 

y  =  1 . 75228763*exp (-alt/6705.6) ; 

else 

y  =  1 . 22557*exp ( -alt/ 9144 ) ; 

end 

end 
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o\°  o\° 


switchlimit  (  u  ) 


function  [  y  ]  = 

SWITCHLIMIT 

Governs  the  output  of  the  switch  in  activating  target  turn  maneuvers. 


O, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 


File : 

Name : 
Compiler : 

Date : 

Description 


Inputs : 
Outputs : 
Process : 
Assumptions 
Comments : 


switchlimit .m 

CPT  Daniel  Perh 

MatLab  v7. 11. 0.584  (R2010b) 

32-bit  (win  32) 

08  July  2011 

In  certain  scenarios,  after  Time-to-Go  falls  below 
specified  threshold  and  the  switch  activates  to 
start  Target  evasive  maneuvers,  the  Time-to-Go  may 
increase  back  above  the  specified  threshold.  When 
this  happens,  the  switch  no  longer  outputs  the 
target  turn  rate  and  the  target  stops  turning  into 
the  missile  for  evasive  maneuvers  until  Time-to-Go 
falls  below  threshold  again. 

Therefore,  this  function  is  implemented  to  output 
the  turn  rate  value  from  the  time  where  the  Time-to- 
Go  falls  below  the  specified  threshold  for  the  first 
time  and  keeps  outputting  the  turn  rate  value 
regardless  if  the  Time-to-Go  subsequently  increases 
above  threshold  value  until  the  end  of  the 
simulation . 

NOTE:  the  global  variables  SWITCHFLAG  and  TURNFLAG 
must  always  be  reset  to  0  prior  to  calling  for  the 
next  model  simulation. 

Target  Turn  Rate 
Target  Turn  Rate 


% -  define  globals  - 

global  SWITCHFLAG  TURNFLAG 

% -  define  constants  - 

% -  define  input  vector  - 

% -  initialize  variables  - 

% -  functions  - 

if  ( (u==0 )  &  ( SWITCHFLAG==0 ) ) 

y=u; 

else  if  (u~=0  &  ( SWITCHFLAG==0 ) ) 
SWITCHFLAG  =  1; 
y  =  u; 

TURNFLAG  =  u; 

else  if  (u~=0  &  ( SWITCHFLAG==1 ) ) 
y=u; 

else  if  (u==0  & ( SWITCHFLAG==1 ) ) 
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y  =  TURNFLAG; 
end 

end 

end 

end 

end 
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function  [  y  ]  =  tgo  (  u  ) 

%  TGO 

%  Computes  time  to  go  from  Range  and  Range  Rate 


O, 

o 

o. 


o 

o, 

o 

File : 

tgo  .m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Computes  Time-to-Go  (tgo) 

o, 

o 

Inputs : 

range,  range  rate 

0, 

o 

Outputs : 

tgo 

o, 

o 

Process : 

o, 

o 

Assumptions : 

o, 

o 

o. 

Comments : 

o 

o, 

o 


% -  define  globals  - 

% -  define  constants  - 

% -  define  input  vector 

range  =  u  (2 ) ; 
rate  =  u ( 1 ) ; 

% -  initialize  variables 

% -  functions  - 

if  (rate==0) 


y  = 

100; 

else 

y  = 

abs (range/rate) ; 

end 

end 
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function  [  y  ]  =  thebigstop  (  u  ) 

%  THEBIGSTOP 

%  Consolidated  simulation  stop  function 


O, 

o 

o. 


o 

o, 

o 

File : 

thebigstop. m 

o, 

o 

Name : 

CPT  Daniel  Perh 

o, 

o 

Compiler : 

MatLab  v7. 11. 0.584  (R2010b) 

o, 

o 

32-bit  (win  32) 

o, 

o 

Date : 

08  July  2011 

o, 

o 

Description : 

Stops  simulation  under  a  variety  of  conditions 

o, 

o 

Inputs : 

see  below 

o, 

o 

Outputs : 

stop  flag 

o, 

o 

Process : 

o, 

o 

Assumptions : 

o, 

o 

CL 

Comments : 

o 

o, 

o 


% -  define  globals  - 

global  STOPFLAG 

% -  define  constants  - 

% - define  input  vector 

R  =  u  ( 1 )  ; 

Rdot  =  u (2) ; 

Vm  =  u ( 3 )  ; 

Vt  =  u ( 4 ) ; 

G  =  u  ( 5 )  ; 

Ny  =  u (6) ; 

N  z  =  u  ( 7  )  ; 
time  =  u  (  8 ) ; 

% -  initialize  variables 

stop  =  [ ]  ; 
y  =  0; 


-  functions  - 

if  ( (time>2 . 0) & (Vm<Vt) ) 
y  =  HI; 

stop  =  ' V  stop ' ; 

end 

if  ( (time>2 . 0 ) & (Rdot>0 ) ) 
y  =  HI; 

stop  =  'Rdot  stop'; 

end 

if  (R<0. 000001) 
y  =  HI; 

stop  =  ' R  stop ' ; 
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end 


if  ( (ST0PFLAG==1) & (y==lll) ) 

disp(['***  ! , stop, '  ***']) 

end 


end 
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APPENDIX  C.  ADDITIONAL  SIMULATION  SCENARIOS 


A  sample  of  the  single  simulation  run  outputs  are  arranged  in  this  appendix  for 
each  of  the  guidance  law  in  a  noise  included  scenario.  For  each  law,  a  different  target  and 
missile  heading  is  chosen  with  different  simulation  parameters  to  demonstrate  the 
capability  of  the  model  and  the  guidance  laws.  The  following  plots  are  then  used  to 
present  the  results: 

•  Engagement  geometry 

•  LOS  range  to  go  vs.  time 

•  Missile  and  target  speed  vs.  time 

•  Missile  and  target  total  acceleration  vs.  time 

•  Missile  forces  (Fx/Fy)  vs.  time 
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A. 


PN  GUIDANCE  LAW 


The  first  scenario  has  an  initial  range  of  20  km,  direct  tail  chase  at  0°  azimuth  but 
with  the  missile  initially  being  launched  pointing  45°  away  from  the  target.  Co-altitude  at 
6000  m,  with  target  6g  maneuver  towards  the  missile  at  three  seconds  tg0. 


Engagement  Geometry 


Missile/Target  Speed  vs  Time 
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Missile  Accelerations  (g)  vs  Time 


c 

■j=  6  -  -  - 

to 

a> 

8  4  - 
o 
< 

0  2  - 
O) 

TO 

I—  Q  _ I _ I _ I _ I _ Lj _ 

0  10  20  30  40  50  60 

Time  (s) 


4 

x  io  Missile  Forces  Fx/Fy  vs  Time 


107 


B. 


APN  GUIDANCE  LAW 


This  second  scenario  has  an  initial  range  of  20  km,  at  45°  azimuth  but  with  the 
missile  initially  being  launched  pointing  -25°  away  from  the  target.  Co-altitude  at  6000 
m,  with  target  6g  maneuver  towards  the  missile  at  three  seconds  tg0. 


Engagement  Geometry 


Missile/Target  Speed  vs  Time 
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Missile  Accelerations  (g)  vs  Time 
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c. 


DG  GUIDANCE  LAW 


This  third  scenario  has  an  initial  range  of  20  km,  at  110°  azimuth  but  with  the 
missile  initially  being  launched  pointing  80°  away  from  the  target.  Co-altitude  at  6000  m. 
In  this  scenario,  the  target  is  modeled  as  moving  with  a  constant  2g  amplitude  sinusoidal 
acceleration  with  a  frequency  of  0.1  n  rad/sec  with  a  final  6g  maneuver  towards  the 
missile  at  three  seconds  tg0. 


Engagement  Geometry 


Missile/Target  Speed  vs  Time 
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Force  (N)  Tarnpt  Ar„lfiratinn  ,n*  Missile  Acceleration  (g) 


Missile  Accelerations  (g)  vs  Time 


Target  Accelerations  (g)  vs  Time 


x  iq4  Missile  Forces  Fx/Fy  vs  Time 
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APPENDIX  D.  NOISE  SIMULATION  RESULTS  DATA 
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Guidance 

Law 

Heading 

Target 

Maneuver 

Max  Range 

Factor  of  Baseline 

Noise 

Simulation  Run  Time 

APN 

45  deg 
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Guidance 

Law 

Heading 

Target 

Maneuver 

Max  Range 

Factor  of  Baseline 

Noise 

Simulation  Run  Time 

DG 

45  deg 

0  turn 

42500  m 

16. 2x 

5214  s 

Min 

0.0002 

Max 

237.1886 

90 

Mean 
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Guidance 

Law 

Heading 

Target 

Maneuver 

Max  Range 

Factor  of  Baseline 

Noise 

Simulation  Run  Time 

PN 

45  deg 

6  turn 

35270  m 

16x 

2797  s 

Min 

0.0019 

Max 

11.3017 

Mean 

3.3956 

Median 

3.1703 

Standard 

Deviation 

2.6143 

Variance 
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Miss 

Distance 
w/o  noise 
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(same  sim 

conditions) 

94  Hits 
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conditions) 

%  Hits 
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Guidance 

Lav/ 

Heading 

Target 

Maneuver 
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Factor  of  Baseline 

Noise 

Simulation  Run  Time 

APN 

45  deg 

6  turn 

42010  m 

7.3x 

4244  s 
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Mean 

Median 

Standard 

Deviation 

Variance 

Miss 
Distance 
w/o  noise 
(same  sim 
conditions) 

94  Hits 
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Guidance 

Law 

Heading 

Target 

Maneuver 

Max  Range 

Factor  of  Baseline 

Noise 

Simulation  Run  Time 

DG 

45  deg 

6  turn 

47030  m 

16. 8x 

4886  s 

Min 

Max 

Mean 

Median 

Standard 

Deviation 

Variance 

Miss 
Distance 
w/o  noise 
(same  sim 
conditions) 

%  Hits 
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Guidance 

Law 

Heading 

Target 

Maneuver 

Max  Range 

Factor  of  Baseline 

Noise 

Simulation  Run  Time 
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Guidance 

Law 

Heading 

Target 

Maneuver 
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Noise 
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